diff --git a/base/comm/psb_cspgather.F90 b/base/comm/psb_cspgather.F90 index e9d3780f3..8a4fe680f 100644 --- a/base/comm/psb_cspgather.F90 +++ b/base/comm/psb_cspgather.F90 @@ -51,13 +51,14 @@ subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_c_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_csp_allgather + + +subroutine psb_lcsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_cspmat_type), intent(inout) :: loca + type(psb_lcspmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_lc_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_c_spk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_c_spk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_lcsp_allgather diff --git a/base/comm/psb_dspgather.F90 b/base/comm/psb_dspgather.F90 index 61416cb1e..83024e104 100644 --- a/base/comm/psb_dspgather.F90 +++ b/base/comm/psb_dspgather.F90 @@ -51,13 +51,14 @@ subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_d_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_dsp_allgather + + +subroutine psb_ldsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_dspmat_type), intent(inout) :: loca + type(psb_ldspmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_ld_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_r_dpk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_r_dpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_ldsp_allgather diff --git a/base/comm/psb_ispgather.F90 b/base/comm/psb_ispgather.F90 index 08d4eea54..3de657b6b 100644 --- a/base/comm/psb_ispgather.F90 +++ b/base/comm/psb_ispgather.F90 @@ -51,13 +51,14 @@ subroutine psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_i_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_isp_allgather + + +subroutine psb_@LX@sp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_ispmat_type), intent(inout) :: loca + type(psb_@LX@spmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_@LX@_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_ipk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_ipk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_@LX@sp_allgather diff --git a/base/comm/psb_lspgather.F90 b/base/comm/psb_lspgather.F90 index e3c1eb051..742c89dea 100644 --- a/base/comm/psb_lspgather.F90 +++ b/base/comm/psb_lspgather.F90 @@ -51,13 +51,14 @@ subroutine psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_l_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_lsp_allgather + + +subroutine psb_@LX@sp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_lspmat_type), intent(inout) :: loca + type(psb_@LX@spmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_@LX@_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_lpk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_@LX@sp_allgather diff --git a/base/comm/psb_sspgather.F90 b/base/comm/psb_sspgather.F90 index 055ebb049..d101b3252 100644 --- a/base/comm/psb_sspgather.F90 +++ b/base/comm/psb_sspgather.F90 @@ -51,13 +51,14 @@ subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_s_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_ssp_allgather + + +subroutine psb_lssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_sspmat_type), intent(inout) :: loca + type(psb_lsspmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_ls_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_r_spk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_r_spk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_lssp_allgather diff --git a/base/comm/psb_zspgather.F90 b/base/comm/psb_zspgather.F90 index 5e62017b9..ba01c9a2b 100644 --- a/base/comm/psb_zspgather.F90 +++ b/base/comm/psb_zspgather.F90 @@ -51,13 +51,14 @@ subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep logical, intent(in), optional :: keepnum,keeploc type(psb_z_coo_sparse_mat) :: loc_coo, glob_coo - integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) - integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg - integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + integer(psb_ipk_) :: nrg, ncg, nzg, nzl + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k logical :: keepnum_, keeploc_ integer(psb_mpk_) :: ictxt,np,me integer(psb_mpk_) :: icomm, minfo, ndx integer(psb_mpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:) integer(psb_ipk_) :: ierr(5) character(len=20) :: name integer(psb_ipk_) :: debug_level, debug_unit @@ -169,3 +170,136 @@ subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep return end subroutine psb_zsp_allgather + + +subroutine psb_lzsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + use psb_desc_mod + use psb_error_mod + use psb_penv_mod + use psb_mat_mod + use psb_tools_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + type(psb_zspmat_type), intent(inout) :: loca + type(psb_lzspmat_type), intent(inout) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root, dupl + logical, intent(in), optional :: keepnum,keeploc + + type(psb_lz_coo_sparse_mat) :: loc_coo, glob_coo + integer(psb_lpk_) :: nrg, ncg, nzg + integer(psb_ipk_) :: err_act, dupl_ + integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl + logical :: keepnum_, keeploc_ + integer(psb_mpk_) :: ictxt,np,me + integer(psb_mpk_) :: icomm, minfo, ndx + integer(psb_lpk_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + integer(psb_ipk_) :: debug_level, debug_unit + + name='psb_gather' + if (psb_get_errstatus().ne.0) return + info=psb_success_ + + call psb_erractionsave(err_act) + ictxt = desc_a%get_context() + icomm = desc_a%get_mpic() + call psb_info(ictxt, me, np) + + if (present(keepnum)) then + keepnum_ = keepnum + else + keepnum_ = .true. + end if + if (present(keeploc)) then + keeploc_ = keeploc + else + keeploc_ = .true. + end if + call globa%free() + + if (keepnum_) then + nrg = desc_a%get_global_rows() + ncg = desc_a%get_global_rows() + + allocate(nzbr(np), idisp(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + ierr(1) = 2*np + call psb_errpush(info,name,i_err=ierr,a_err='integer') + goto 9999 + end if + + + if (keeploc_) then + call loca%cp_to(loc_coo) + else + call loca%mv_to(loc_coo) + end if + nzl = loc_coo%get_nzeros() + call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') + call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') + nzbr(:) = 0 + nzbr(me+1) = nzl + call psb_sum(ictxt,nzbr(1:np)) + nzg = sum(nzbr) + if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg) + if (info /= psb_success_) goto 9999 + + do ip=1,np + idisp(ip) = sum(nzbr(1:ip-1)) + enddo + ndx = nzbr(me+1) + call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_c_dpk_,& + & glob_coo%val,nzbr,idisp,& + & psb_mpi_c_dpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,& + & glob_coo%ia,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + if (minfo == psb_success_) call & + & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,& + & glob_coo%ja,nzbr,idisp,& + & psb_mpi_lpk_,icomm,minfo) + + if (minfo /= psb_success_) then + info = minfo + call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv') + goto 9999 + end if + call loc_coo%free() + ! + ! Is the code below safe? For very large cases + ! the indices in glob_coo will overflow. But then, + ! for very large cases it does not make sense to + ! gather the matrix on a single procecss anyway... + ! + call glob_coo%set_nzeros(nzg) + if (present(dupl)) call glob_coo%set_dupl(dupl) + call globa%mv_from(glob_coo) + + else + write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_ + info = -1 + goto 9999 + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_error_handler(ione*ictxt,err_act) + + return + +end subroutine psb_lzsp_allgather diff --git a/base/modules/comm/psb_c_comm_mod.f90 b/base/modules/comm/psb_c_comm_mod.f90 index 544e8de5a..0deb379f9 100644 --- a/base/modules/comm/psb_c_comm_mod.f90 +++ b/base/modules/comm/psb_c_comm_mod.f90 @@ -31,7 +31,7 @@ ! module psb_c_comm_mod use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_spk_ - use psb_mat_mod, only : psb_cspmat_type + use psb_mat_mod, only : psb_cspmat_type, psb_lcspmat_type use psb_c_vect_mod, only : psb_c_vect_type, psb_c_base_vect_type use psb_c_multivect_mod, only : psb_c_multivect_type, psb_c_base_multivect_type @@ -105,6 +105,16 @@ module psb_c_comm_mod integer(psb_ipk_), intent(in), optional :: root,dupl logical, intent(in), optional :: keepnum,keeploc end subroutine psb_csp_allgather + subroutine psb_lcsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + import + implicit none + type(psb_cspmat_type), intent(inout) :: loca + type(psb_lcspmat_type), intent(out) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root,dupl + logical, intent(in), optional :: keepnum,keeploc + end subroutine psb_lcsp_allgather subroutine psb_cgather_vect(globx, locx, desc_a, info, root) import implicit none diff --git a/base/modules/comm/psb_d_comm_mod.f90 b/base/modules/comm/psb_d_comm_mod.f90 index 2b4bf4b0d..431a8aec0 100644 --- a/base/modules/comm/psb_d_comm_mod.f90 +++ b/base/modules/comm/psb_d_comm_mod.f90 @@ -31,7 +31,7 @@ ! module psb_d_comm_mod use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_dpk_ - use psb_mat_mod, only : psb_dspmat_type + use psb_mat_mod, only : psb_dspmat_type, psb_ldspmat_type use psb_d_vect_mod, only : psb_d_vect_type, psb_d_base_vect_type use psb_d_multivect_mod, only : psb_d_multivect_type, psb_d_base_multivect_type @@ -105,6 +105,16 @@ module psb_d_comm_mod integer(psb_ipk_), intent(in), optional :: root,dupl logical, intent(in), optional :: keepnum,keeploc end subroutine psb_dsp_allgather + subroutine psb_ldsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + import + implicit none + type(psb_dspmat_type), intent(inout) :: loca + type(psb_ldspmat_type), intent(out) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root,dupl + logical, intent(in), optional :: keepnum,keeploc + end subroutine psb_ldsp_allgather subroutine psb_dgather_vect(globx, locx, desc_a, info, root) import implicit none diff --git a/base/modules/comm/psb_s_comm_mod.f90 b/base/modules/comm/psb_s_comm_mod.f90 index 77fa92d91..5e54e360b 100644 --- a/base/modules/comm/psb_s_comm_mod.f90 +++ b/base/modules/comm/psb_s_comm_mod.f90 @@ -31,7 +31,7 @@ ! module psb_s_comm_mod use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_spk_ - use psb_mat_mod, only : psb_sspmat_type + use psb_mat_mod, only : psb_sspmat_type, psb_lsspmat_type use psb_s_vect_mod, only : psb_s_vect_type, psb_s_base_vect_type use psb_s_multivect_mod, only : psb_s_multivect_type, psb_s_base_multivect_type @@ -105,6 +105,16 @@ module psb_s_comm_mod integer(psb_ipk_), intent(in), optional :: root,dupl logical, intent(in), optional :: keepnum,keeploc end subroutine psb_ssp_allgather + subroutine psb_lssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + import + implicit none + type(psb_sspmat_type), intent(inout) :: loca + type(psb_lsspmat_type), intent(out) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root,dupl + logical, intent(in), optional :: keepnum,keeploc + end subroutine psb_lssp_allgather subroutine psb_sgather_vect(globx, locx, desc_a, info, root) import implicit none diff --git a/base/modules/comm/psb_z_comm_mod.f90 b/base/modules/comm/psb_z_comm_mod.f90 index 5367ec10d..bb38d7253 100644 --- a/base/modules/comm/psb_z_comm_mod.f90 +++ b/base/modules/comm/psb_z_comm_mod.f90 @@ -31,7 +31,7 @@ ! module psb_z_comm_mod use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_dpk_ - use psb_mat_mod, only : psb_zspmat_type + use psb_mat_mod, only : psb_zspmat_type, psb_lzspmat_type use psb_z_vect_mod, only : psb_z_vect_type, psb_z_base_vect_type use psb_z_multivect_mod, only : psb_z_multivect_type, psb_z_base_multivect_type @@ -105,6 +105,16 @@ module psb_z_comm_mod integer(psb_ipk_), intent(in), optional :: root,dupl logical, intent(in), optional :: keepnum,keeploc end subroutine psb_zsp_allgather + subroutine psb_lzsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) + import + implicit none + type(psb_zspmat_type), intent(inout) :: loca + type(psb_lzspmat_type), intent(out) :: globa + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: root,dupl + logical, intent(in), optional :: keepnum,keeploc + end subroutine psb_lzsp_allgather subroutine psb_zgather_vect(globx, locx, desc_a, info, root) import implicit none diff --git a/base/modules/serial/psb_c_mat_mod.F90 b/base/modules/serial/psb_c_mat_mod.F90 index 5756df426..442f1322e 100644 --- a/base/modules/serial/psb_c_mat_mod.F90 +++ b/base/modules/serial/psb_c_mat_mod.F90 @@ -194,7 +194,22 @@ module psb_c_mat_mod procedure, pass(a) :: cscnv_base => psb_c_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_cspmat_clone - + ! + ! To/from lc + ! + procedure, pass(a) :: mv_from_lb => psb_c_mv_from_lb + procedure, pass(a) :: mv_to_lb => psb_c_mv_to_lb + procedure, pass(a) :: cp_from_lb => psb_c_cp_from_lb + procedure, pass(a) :: cp_to_lb => psb_c_cp_to_lb + procedure, pass(a) :: mv_from_l => psb_c_mv_from_l + procedure, pass(a) :: mv_to_l => psb_c_mv_to_l + procedure, pass(a) :: cp_from_l => psb_c_cp_from_l + procedure, pass(a) :: cp_to_l => psb_c_cp_to_l + generic, public :: mv_from => mv_from_lb, mv_from_l + generic, public :: mv_to => mv_to_lb, mv_to_l + generic, public :: cp_from => cp_from_lb, cp_from_l + generic, public :: cp_to => cp_to_lb, cp_to_l + ! Computational routines procedure, pass(a) :: get_diag => psb_c_get_diag procedure, pass(a) :: maxval => psb_c_maxval @@ -355,6 +370,21 @@ module psb_c_mat_mod procedure, pass(a) :: cscnv_base => psb_lc_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_lcspmat_clone + ! + ! To/from c + ! + procedure, pass(a) :: mv_from_ib => psb_lc_mv_from_ib + procedure, pass(a) :: mv_to_ib => psb_lc_mv_to_ib + procedure, pass(a) :: cp_from_ib => psb_lc_cp_from_ib + procedure, pass(a) :: cp_to_ib => psb_lc_cp_to_ib + procedure, pass(a) :: mv_from_i => psb_lc_mv_from_i + procedure, pass(a) :: mv_to_i => psb_lc_mv_to_i + procedure, pass(a) :: cp_from_i => psb_lc_cp_from_i + procedure, pass(a) :: cp_to_i => psb_lc_cp_to_i + generic, public :: mv_from => mv_from_ib, mv_from_i + generic, public :: mv_to => mv_to_ib, mv_to_i + generic, public :: cp_from => cp_from_ib, cp_from_i + generic, public :: cp_to => cp_to_ib, cp_to_i ! Computational routines procedure, pass(a) :: get_diag => psb_lc_get_diag @@ -845,7 +875,73 @@ module psb_c_mat_mod class(psb_c_base_sparse_mat), intent(inout) :: b end subroutine psb_c_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_c_mv_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat + class(psb_cspmat_type), intent(inout) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + end subroutine psb_c_mv_from_lb + end interface + + interface + subroutine psb_c_cp_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat + class(psb_cspmat_type), intent(out) :: a + class(psb_lc_base_sparse_mat), intent(in) :: b + end subroutine psb_c_cp_from_lb + end interface + + interface + subroutine psb_c_mv_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat + class(psb_cspmat_type), intent(inout) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + end subroutine psb_c_mv_to_lb + end interface + interface + subroutine psb_c_cp_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat + class(psb_cspmat_type), intent(in) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + end subroutine psb_c_cp_to_lb + end interface + + interface + subroutine psb_c_mv_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type + class(psb_cspmat_type), intent(inout) :: a + class(psb_lcspmat_type), intent(inout) :: b + end subroutine psb_c_mv_from_l + end interface + + interface + subroutine psb_c_cp_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type + class(psb_cspmat_type), intent(out) :: a + class(psb_lcspmat_type), intent(in) :: b + end subroutine psb_c_cp_from_l + end interface + + interface + subroutine psb_c_mv_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type + class(psb_cspmat_type), intent(inout) :: a + class(psb_lcspmat_type), intent(inout) :: b + end subroutine psb_c_mv_to_l + end interface + + interface + subroutine psb_c_cp_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type + class(psb_cspmat_type), intent(in) :: a + class(psb_lcspmat_type), intent(inout) :: b + end subroutine psb_c_cp_to_l + end interface + ! ! Transfer the internal allocation to the target. ! @@ -1467,6 +1563,73 @@ module psb_c_mat_mod class(psb_lc_base_sparse_mat), intent(inout) :: b end subroutine psb_lc_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_lc_mv_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat + class(psb_lcspmat_type), intent(inout) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + end subroutine psb_lc_mv_from_ib + end interface + + interface + subroutine psb_lc_cp_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat + class(psb_lcspmat_type), intent(out) :: a + class(psb_c_base_sparse_mat), intent(in) :: b + end subroutine psb_lc_cp_from_ib + end interface + + interface + subroutine psb_lc_mv_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat + class(psb_lcspmat_type), intent(inout) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + end subroutine psb_lc_mv_to_ib + end interface + + interface + subroutine psb_lc_cp_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat + class(psb_lcspmat_type), intent(in) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + end subroutine psb_lc_cp_to_ib + end interface + + interface + subroutine psb_lc_mv_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type + class(psb_lcspmat_type), intent(inout) :: a + class(psb_cspmat_type), intent(inout) :: b + end subroutine psb_lc_mv_from_i + end interface + + interface + subroutine psb_lc_cp_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type + class(psb_lcspmat_type), intent(out) :: a + class(psb_cspmat_type), intent(in) :: b + end subroutine psb_lc_cp_from_i + end interface + + interface + subroutine psb_lc_mv_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type + class(psb_lcspmat_type), intent(inout) :: a + class(psb_cspmat_type), intent(inout) :: b + end subroutine psb_lc_mv_to_i + end interface + + interface + subroutine psb_lc_cp_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type + class(psb_lcspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + end subroutine psb_lc_cp_to_i + end interface + ! ! Transfer the internal allocation to the target. diff --git a/base/modules/serial/psb_c_serial_mod.f90 b/base/modules/serial/psb_c_serial_mod.f90 index 5e25a0354..b385e03bd 100644 --- a/base/modules/serial/psb_c_serial_mod.f90 +++ b/base/modules/serial/psb_c_serial_mod.f90 @@ -204,14 +204,14 @@ module psb_c_serial_mod end interface psb_aspxpby interface psb_spspmm -!!$ subroutine psb_cspspmm(a,b,c,info) -!!$ use psb_c_mat_mod, only : psb_cspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_cspmat_type), intent(in) :: a,b -!!$ type(psb_cspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_cspspmm + subroutine psb_lcspspmm(a,b,c,info) + use psb_c_mat_mod, only : psb_lcspmat_type + import :: psb_ipk_ + implicit none + type(psb_lcspmat_type), intent(in) :: a,b + type(psb_lcspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lcspspmm subroutine psb_lccsrspspmm(a,b,c,info) use psb_c_mat_mod, only : psb_lc_csr_sparse_mat import :: psb_ipk_ @@ -231,14 +231,14 @@ module psb_c_serial_mod end interface psb_spspmm interface psb_symbmm -!!$ subroutine psb_csymbmm(a,b,c,info) -!!$ use psb_c_mat_mod, only : psb_cspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_cspmat_type), intent(in) :: a,b -!!$ type(psb_cspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_csymbmm + subroutine psb_lcsymbmm(a,b,c,info) + use psb_c_mat_mod, only : psb_lcspmat_type + import :: psb_ipk_ + implicit none + type(psb_lcspmat_type), intent(in) :: a,b + type(psb_lcspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lcsymbmm subroutine psb_lcbase_symbmm(a,b,c,info) use psb_c_mat_mod, only : psb_lc_base_sparse_mat, psb_lc_csr_sparse_mat import :: psb_ipk_ @@ -250,13 +250,13 @@ module psb_c_serial_mod end interface psb_symbmm interface psb_numbmm -!!$ subroutine psb_cnumbmm(a,b,c) -!!$ use psb_c_mat_mod, only : psb_cspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_cspmat_type), intent(in) :: a,b -!!$ type(psb_cspmat_type), intent(inout) :: c -!!$ end subroutine psb_cnumbmm + subroutine psb_lcnumbmm(a,b,c) + use psb_c_mat_mod, only : psb_lcspmat_type + import :: psb_ipk_ + implicit none + type(psb_lcspmat_type), intent(in) :: a,b + type(psb_lcspmat_type), intent(inout) :: c + end subroutine psb_lcnumbmm subroutine psb_lcbase_numbmm(a,b,c) use psb_c_mat_mod, only : psb_lc_base_sparse_mat, psb_lc_csr_sparse_mat import :: psb_ipk_ diff --git a/base/modules/serial/psb_d_mat_mod.F90 b/base/modules/serial/psb_d_mat_mod.F90 index 76b5eae83..acc75baf4 100644 --- a/base/modules/serial/psb_d_mat_mod.F90 +++ b/base/modules/serial/psb_d_mat_mod.F90 @@ -194,7 +194,22 @@ module psb_d_mat_mod procedure, pass(a) :: cscnv_base => psb_d_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_dspmat_clone - + ! + ! To/from ld + ! + procedure, pass(a) :: mv_from_lb => psb_d_mv_from_lb + procedure, pass(a) :: mv_to_lb => psb_d_mv_to_lb + procedure, pass(a) :: cp_from_lb => psb_d_cp_from_lb + procedure, pass(a) :: cp_to_lb => psb_d_cp_to_lb + procedure, pass(a) :: mv_from_l => psb_d_mv_from_l + procedure, pass(a) :: mv_to_l => psb_d_mv_to_l + procedure, pass(a) :: cp_from_l => psb_d_cp_from_l + procedure, pass(a) :: cp_to_l => psb_d_cp_to_l + generic, public :: mv_from => mv_from_lb, mv_from_l + generic, public :: mv_to => mv_to_lb, mv_to_l + generic, public :: cp_from => cp_from_lb, cp_from_l + generic, public :: cp_to => cp_to_lb, cp_to_l + ! Computational routines procedure, pass(a) :: get_diag => psb_d_get_diag procedure, pass(a) :: maxval => psb_d_maxval @@ -355,6 +370,21 @@ module psb_d_mat_mod procedure, pass(a) :: cscnv_base => psb_ld_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_ldspmat_clone + ! + ! To/from d + ! + procedure, pass(a) :: mv_from_ib => psb_ld_mv_from_ib + procedure, pass(a) :: mv_to_ib => psb_ld_mv_to_ib + procedure, pass(a) :: cp_from_ib => psb_ld_cp_from_ib + procedure, pass(a) :: cp_to_ib => psb_ld_cp_to_ib + procedure, pass(a) :: mv_from_i => psb_ld_mv_from_i + procedure, pass(a) :: mv_to_i => psb_ld_mv_to_i + procedure, pass(a) :: cp_from_i => psb_ld_cp_from_i + procedure, pass(a) :: cp_to_i => psb_ld_cp_to_i + generic, public :: mv_from => mv_from_ib, mv_from_i + generic, public :: mv_to => mv_to_ib, mv_to_i + generic, public :: cp_from => cp_from_ib, cp_from_i + generic, public :: cp_to => cp_to_ib, cp_to_i ! Computational routines procedure, pass(a) :: get_diag => psb_ld_get_diag @@ -845,7 +875,73 @@ module psb_d_mat_mod class(psb_d_base_sparse_mat), intent(inout) :: b end subroutine psb_d_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_d_mv_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat + class(psb_dspmat_type), intent(inout) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + end subroutine psb_d_mv_from_lb + end interface + + interface + subroutine psb_d_cp_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat + class(psb_dspmat_type), intent(out) :: a + class(psb_ld_base_sparse_mat), intent(in) :: b + end subroutine psb_d_cp_from_lb + end interface + + interface + subroutine psb_d_mv_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat + class(psb_dspmat_type), intent(inout) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + end subroutine psb_d_mv_to_lb + end interface + interface + subroutine psb_d_cp_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat + class(psb_dspmat_type), intent(in) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + end subroutine psb_d_cp_to_lb + end interface + + interface + subroutine psb_d_mv_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type + class(psb_dspmat_type), intent(inout) :: a + class(psb_ldspmat_type), intent(inout) :: b + end subroutine psb_d_mv_from_l + end interface + + interface + subroutine psb_d_cp_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type + class(psb_dspmat_type), intent(out) :: a + class(psb_ldspmat_type), intent(in) :: b + end subroutine psb_d_cp_from_l + end interface + + interface + subroutine psb_d_mv_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type + class(psb_dspmat_type), intent(inout) :: a + class(psb_ldspmat_type), intent(inout) :: b + end subroutine psb_d_mv_to_l + end interface + + interface + subroutine psb_d_cp_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type + class(psb_dspmat_type), intent(in) :: a + class(psb_ldspmat_type), intent(inout) :: b + end subroutine psb_d_cp_to_l + end interface + ! ! Transfer the internal allocation to the target. ! @@ -1467,6 +1563,73 @@ module psb_d_mat_mod class(psb_ld_base_sparse_mat), intent(inout) :: b end subroutine psb_ld_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_ld_mv_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat + class(psb_ldspmat_type), intent(inout) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + end subroutine psb_ld_mv_from_ib + end interface + + interface + subroutine psb_ld_cp_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat + class(psb_ldspmat_type), intent(out) :: a + class(psb_d_base_sparse_mat), intent(in) :: b + end subroutine psb_ld_cp_from_ib + end interface + + interface + subroutine psb_ld_mv_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat + class(psb_ldspmat_type), intent(inout) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + end subroutine psb_ld_mv_to_ib + end interface + + interface + subroutine psb_ld_cp_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat + class(psb_ldspmat_type), intent(in) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + end subroutine psb_ld_cp_to_ib + end interface + + interface + subroutine psb_ld_mv_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type + class(psb_ldspmat_type), intent(inout) :: a + class(psb_dspmat_type), intent(inout) :: b + end subroutine psb_ld_mv_from_i + end interface + + interface + subroutine psb_ld_cp_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type + class(psb_ldspmat_type), intent(out) :: a + class(psb_dspmat_type), intent(in) :: b + end subroutine psb_ld_cp_from_i + end interface + + interface + subroutine psb_ld_mv_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type + class(psb_ldspmat_type), intent(inout) :: a + class(psb_dspmat_type), intent(inout) :: b + end subroutine psb_ld_mv_to_i + end interface + + interface + subroutine psb_ld_cp_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type + class(psb_ldspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + end subroutine psb_ld_cp_to_i + end interface + ! ! Transfer the internal allocation to the target. diff --git a/base/modules/serial/psb_d_serial_mod.f90 b/base/modules/serial/psb_d_serial_mod.f90 index 84b4781f0..7c45d30d2 100644 --- a/base/modules/serial/psb_d_serial_mod.f90 +++ b/base/modules/serial/psb_d_serial_mod.f90 @@ -204,14 +204,14 @@ module psb_d_serial_mod end interface psb_aspxpby interface psb_spspmm -!!$ subroutine psb_dspspmm(a,b,c,info) -!!$ use psb_d_mat_mod, only : psb_dspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_dspmat_type), intent(in) :: a,b -!!$ type(psb_dspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_dspspmm + subroutine psb_ldspspmm(a,b,c,info) + use psb_d_mat_mod, only : psb_ldspmat_type + import :: psb_ipk_ + implicit none + type(psb_ldspmat_type), intent(in) :: a,b + type(psb_ldspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_ldspspmm subroutine psb_ldcsrspspmm(a,b,c,info) use psb_d_mat_mod, only : psb_ld_csr_sparse_mat import :: psb_ipk_ @@ -231,14 +231,14 @@ module psb_d_serial_mod end interface psb_spspmm interface psb_symbmm -!!$ subroutine psb_dsymbmm(a,b,c,info) -!!$ use psb_d_mat_mod, only : psb_dspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_dspmat_type), intent(in) :: a,b -!!$ type(psb_dspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_dsymbmm + subroutine psb_ldsymbmm(a,b,c,info) + use psb_d_mat_mod, only : psb_ldspmat_type + import :: psb_ipk_ + implicit none + type(psb_ldspmat_type), intent(in) :: a,b + type(psb_ldspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_ldsymbmm subroutine psb_ldbase_symbmm(a,b,c,info) use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat import :: psb_ipk_ @@ -250,13 +250,13 @@ module psb_d_serial_mod end interface psb_symbmm interface psb_numbmm -!!$ subroutine psb_dnumbmm(a,b,c) -!!$ use psb_d_mat_mod, only : psb_dspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_dspmat_type), intent(in) :: a,b -!!$ type(psb_dspmat_type), intent(inout) :: c -!!$ end subroutine psb_dnumbmm + subroutine psb_ldnumbmm(a,b,c) + use psb_d_mat_mod, only : psb_ldspmat_type + import :: psb_ipk_ + implicit none + type(psb_ldspmat_type), intent(in) :: a,b + type(psb_ldspmat_type), intent(inout) :: c + end subroutine psb_ldnumbmm subroutine psb_ldbase_numbmm(a,b,c) use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat import :: psb_ipk_ diff --git a/base/modules/serial/psb_s_mat_mod.F90 b/base/modules/serial/psb_s_mat_mod.F90 index c8d036960..b1b5eb24b 100644 --- a/base/modules/serial/psb_s_mat_mod.F90 +++ b/base/modules/serial/psb_s_mat_mod.F90 @@ -194,7 +194,22 @@ module psb_s_mat_mod procedure, pass(a) :: cscnv_base => psb_s_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_sspmat_clone - + ! + ! To/from ls + ! + procedure, pass(a) :: mv_from_lb => psb_s_mv_from_lb + procedure, pass(a) :: mv_to_lb => psb_s_mv_to_lb + procedure, pass(a) :: cp_from_lb => psb_s_cp_from_lb + procedure, pass(a) :: cp_to_lb => psb_s_cp_to_lb + procedure, pass(a) :: mv_from_l => psb_s_mv_from_l + procedure, pass(a) :: mv_to_l => psb_s_mv_to_l + procedure, pass(a) :: cp_from_l => psb_s_cp_from_l + procedure, pass(a) :: cp_to_l => psb_s_cp_to_l + generic, public :: mv_from => mv_from_lb, mv_from_l + generic, public :: mv_to => mv_to_lb, mv_to_l + generic, public :: cp_from => cp_from_lb, cp_from_l + generic, public :: cp_to => cp_to_lb, cp_to_l + ! Computational routines procedure, pass(a) :: get_diag => psb_s_get_diag procedure, pass(a) :: maxval => psb_s_maxval @@ -355,6 +370,21 @@ module psb_s_mat_mod procedure, pass(a) :: cscnv_base => psb_ls_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_lsspmat_clone + ! + ! To/from s + ! + procedure, pass(a) :: mv_from_ib => psb_ls_mv_from_ib + procedure, pass(a) :: mv_to_ib => psb_ls_mv_to_ib + procedure, pass(a) :: cp_from_ib => psb_ls_cp_from_ib + procedure, pass(a) :: cp_to_ib => psb_ls_cp_to_ib + procedure, pass(a) :: mv_from_i => psb_ls_mv_from_i + procedure, pass(a) :: mv_to_i => psb_ls_mv_to_i + procedure, pass(a) :: cp_from_i => psb_ls_cp_from_i + procedure, pass(a) :: cp_to_i => psb_ls_cp_to_i + generic, public :: mv_from => mv_from_ib, mv_from_i + generic, public :: mv_to => mv_to_ib, mv_to_i + generic, public :: cp_from => cp_from_ib, cp_from_i + generic, public :: cp_to => cp_to_ib, cp_to_i ! Computational routines procedure, pass(a) :: get_diag => psb_ls_get_diag @@ -845,7 +875,73 @@ module psb_s_mat_mod class(psb_s_base_sparse_mat), intent(inout) :: b end subroutine psb_s_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_s_mv_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat + class(psb_sspmat_type), intent(inout) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + end subroutine psb_s_mv_from_lb + end interface + + interface + subroutine psb_s_cp_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat + class(psb_sspmat_type), intent(out) :: a + class(psb_ls_base_sparse_mat), intent(in) :: b + end subroutine psb_s_cp_from_lb + end interface + + interface + subroutine psb_s_mv_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat + class(psb_sspmat_type), intent(inout) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + end subroutine psb_s_mv_to_lb + end interface + interface + subroutine psb_s_cp_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat + class(psb_sspmat_type), intent(in) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + end subroutine psb_s_cp_to_lb + end interface + + interface + subroutine psb_s_mv_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type + class(psb_sspmat_type), intent(inout) :: a + class(psb_lsspmat_type), intent(inout) :: b + end subroutine psb_s_mv_from_l + end interface + + interface + subroutine psb_s_cp_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type + class(psb_sspmat_type), intent(out) :: a + class(psb_lsspmat_type), intent(in) :: b + end subroutine psb_s_cp_from_l + end interface + + interface + subroutine psb_s_mv_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type + class(psb_sspmat_type), intent(inout) :: a + class(psb_lsspmat_type), intent(inout) :: b + end subroutine psb_s_mv_to_l + end interface + + interface + subroutine psb_s_cp_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type + class(psb_sspmat_type), intent(in) :: a + class(psb_lsspmat_type), intent(inout) :: b + end subroutine psb_s_cp_to_l + end interface + ! ! Transfer the internal allocation to the target. ! @@ -1467,6 +1563,73 @@ module psb_s_mat_mod class(psb_ls_base_sparse_mat), intent(inout) :: b end subroutine psb_ls_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_ls_mv_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat + class(psb_lsspmat_type), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + end subroutine psb_ls_mv_from_ib + end interface + + interface + subroutine psb_ls_cp_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat + class(psb_lsspmat_type), intent(out) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + end subroutine psb_ls_cp_from_ib + end interface + + interface + subroutine psb_ls_mv_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat + class(psb_lsspmat_type), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + end subroutine psb_ls_mv_to_ib + end interface + + interface + subroutine psb_ls_cp_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat + class(psb_lsspmat_type), intent(in) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + end subroutine psb_ls_cp_to_ib + end interface + + interface + subroutine psb_ls_mv_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type + class(psb_lsspmat_type), intent(inout) :: a + class(psb_sspmat_type), intent(inout) :: b + end subroutine psb_ls_mv_from_i + end interface + + interface + subroutine psb_ls_cp_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type + class(psb_lsspmat_type), intent(out) :: a + class(psb_sspmat_type), intent(in) :: b + end subroutine psb_ls_cp_from_i + end interface + + interface + subroutine psb_ls_mv_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type + class(psb_lsspmat_type), intent(inout) :: a + class(psb_sspmat_type), intent(inout) :: b + end subroutine psb_ls_mv_to_i + end interface + + interface + subroutine psb_ls_cp_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type + class(psb_lsspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + end subroutine psb_ls_cp_to_i + end interface + ! ! Transfer the internal allocation to the target. diff --git a/base/modules/serial/psb_s_serial_mod.f90 b/base/modules/serial/psb_s_serial_mod.f90 index 09f4a0f7e..0e872e0df 100644 --- a/base/modules/serial/psb_s_serial_mod.f90 +++ b/base/modules/serial/psb_s_serial_mod.f90 @@ -204,14 +204,14 @@ module psb_s_serial_mod end interface psb_aspxpby interface psb_spspmm -!!$ subroutine psb_sspspmm(a,b,c,info) -!!$ use psb_s_mat_mod, only : psb_sspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_sspmat_type), intent(in) :: a,b -!!$ type(psb_sspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_sspspmm + subroutine psb_lsspspmm(a,b,c,info) + use psb_s_mat_mod, only : psb_lsspmat_type + import :: psb_ipk_ + implicit none + type(psb_lsspmat_type), intent(in) :: a,b + type(psb_lsspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lsspspmm subroutine psb_lscsrspspmm(a,b,c,info) use psb_s_mat_mod, only : psb_ls_csr_sparse_mat import :: psb_ipk_ @@ -231,14 +231,14 @@ module psb_s_serial_mod end interface psb_spspmm interface psb_symbmm -!!$ subroutine psb_ssymbmm(a,b,c,info) -!!$ use psb_s_mat_mod, only : psb_sspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_sspmat_type), intent(in) :: a,b -!!$ type(psb_sspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_ssymbmm + subroutine psb_lssymbmm(a,b,c,info) + use psb_s_mat_mod, only : psb_lsspmat_type + import :: psb_ipk_ + implicit none + type(psb_lsspmat_type), intent(in) :: a,b + type(psb_lsspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lssymbmm subroutine psb_lsbase_symbmm(a,b,c,info) use psb_s_mat_mod, only : psb_ls_base_sparse_mat, psb_ls_csr_sparse_mat import :: psb_ipk_ @@ -250,13 +250,13 @@ module psb_s_serial_mod end interface psb_symbmm interface psb_numbmm -!!$ subroutine psb_snumbmm(a,b,c) -!!$ use psb_s_mat_mod, only : psb_sspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_sspmat_type), intent(in) :: a,b -!!$ type(psb_sspmat_type), intent(inout) :: c -!!$ end subroutine psb_snumbmm + subroutine psb_lsnumbmm(a,b,c) + use psb_s_mat_mod, only : psb_lsspmat_type + import :: psb_ipk_ + implicit none + type(psb_lsspmat_type), intent(in) :: a,b + type(psb_lsspmat_type), intent(inout) :: c + end subroutine psb_lsnumbmm subroutine psb_lsbase_numbmm(a,b,c) use psb_s_mat_mod, only : psb_ls_base_sparse_mat, psb_ls_csr_sparse_mat import :: psb_ipk_ diff --git a/base/modules/serial/psb_serial_mod.f90 b/base/modules/serial/psb_serial_mod.f90 index 97bff01e8..2f2154e0e 100644 --- a/base/modules/serial/psb_serial_mod.f90 +++ b/base/modules/serial/psb_serial_mod.f90 @@ -68,7 +68,7 @@ module psb_serial_mod end subroutine psb_d_nspaxpby end interface psb_nspaxpby - interface symbmm + interface subroutine symbmm (n, m, l, ia, ja, diaga, & & ib, jb, diagb, ic, jc, diagc, index) import :: psb_ipk_ @@ -76,6 +76,13 @@ module psb_serial_mod & diagc, index(*) integer(psb_ipk_), allocatable :: ic(:),jc(:) end subroutine symbmm + subroutine lsymbmm (n, m, l, ia, ja, diaga, & + & ib, jb, diagb, ic, jc, diagc, index) + import :: psb_ipk_, psb_lpk_ + integer(psb_lpk_) :: n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& + & diagc, index(*) + integer(psb_lpk_), allocatable :: ic(:),jc(:) + end subroutine lsymbmm end interface diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90 index 9e6b253cd..f40d35cb7 100644 --- a/base/modules/serial/psb_z_mat_mod.F90 +++ b/base/modules/serial/psb_z_mat_mod.F90 @@ -194,7 +194,22 @@ module psb_z_mat_mod procedure, pass(a) :: cscnv_base => psb_z_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_zspmat_clone - + ! + ! To/from lz + ! + procedure, pass(a) :: mv_from_lb => psb_z_mv_from_lb + procedure, pass(a) :: mv_to_lb => psb_z_mv_to_lb + procedure, pass(a) :: cp_from_lb => psb_z_cp_from_lb + procedure, pass(a) :: cp_to_lb => psb_z_cp_to_lb + procedure, pass(a) :: mv_from_l => psb_z_mv_from_l + procedure, pass(a) :: mv_to_l => psb_z_mv_to_l + procedure, pass(a) :: cp_from_l => psb_z_cp_from_l + procedure, pass(a) :: cp_to_l => psb_z_cp_to_l + generic, public :: mv_from => mv_from_lb, mv_from_l + generic, public :: mv_to => mv_to_lb, mv_to_l + generic, public :: cp_from => cp_from_lb, cp_from_l + generic, public :: cp_to => cp_to_lb, cp_to_l + ! Computational routines procedure, pass(a) :: get_diag => psb_z_get_diag procedure, pass(a) :: maxval => psb_z_maxval @@ -355,6 +370,21 @@ module psb_z_mat_mod procedure, pass(a) :: cscnv_base => psb_lz_cscnv_base generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base procedure, pass(a) :: clone => psb_lzspmat_clone + ! + ! To/from z + ! + procedure, pass(a) :: mv_from_ib => psb_lz_mv_from_ib + procedure, pass(a) :: mv_to_ib => psb_lz_mv_to_ib + procedure, pass(a) :: cp_from_ib => psb_lz_cp_from_ib + procedure, pass(a) :: cp_to_ib => psb_lz_cp_to_ib + procedure, pass(a) :: mv_from_i => psb_lz_mv_from_i + procedure, pass(a) :: mv_to_i => psb_lz_mv_to_i + procedure, pass(a) :: cp_from_i => psb_lz_cp_from_i + procedure, pass(a) :: cp_to_i => psb_lz_cp_to_i + generic, public :: mv_from => mv_from_ib, mv_from_i + generic, public :: mv_to => mv_to_ib, mv_to_i + generic, public :: cp_from => cp_from_ib, cp_from_i + generic, public :: cp_to => cp_to_ib, cp_to_i ! Computational routines procedure, pass(a) :: get_diag => psb_lz_get_diag @@ -845,7 +875,73 @@ module psb_z_mat_mod class(psb_z_base_sparse_mat), intent(inout) :: b end subroutine psb_z_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_z_mv_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat + class(psb_zspmat_type), intent(inout) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + end subroutine psb_z_mv_from_lb + end interface + + interface + subroutine psb_z_cp_from_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat + class(psb_zspmat_type), intent(out) :: a + class(psb_lz_base_sparse_mat), intent(in) :: b + end subroutine psb_z_cp_from_lb + end interface + + interface + subroutine psb_z_mv_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat + class(psb_zspmat_type), intent(inout) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + end subroutine psb_z_mv_to_lb + end interface + interface + subroutine psb_z_cp_to_lb(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat + class(psb_zspmat_type), intent(in) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + end subroutine psb_z_cp_to_lb + end interface + + interface + subroutine psb_z_mv_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type + class(psb_zspmat_type), intent(inout) :: a + class(psb_lzspmat_type), intent(inout) :: b + end subroutine psb_z_mv_from_l + end interface + + interface + subroutine psb_z_cp_from_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type + class(psb_zspmat_type), intent(out) :: a + class(psb_lzspmat_type), intent(in) :: b + end subroutine psb_z_cp_from_l + end interface + + interface + subroutine psb_z_mv_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type + class(psb_zspmat_type), intent(inout) :: a + class(psb_lzspmat_type), intent(inout) :: b + end subroutine psb_z_mv_to_l + end interface + + interface + subroutine psb_z_cp_to_l(a,b) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type + class(psb_zspmat_type), intent(in) :: a + class(psb_lzspmat_type), intent(inout) :: b + end subroutine psb_z_cp_to_l + end interface + ! ! Transfer the internal allocation to the target. ! @@ -1467,6 +1563,73 @@ module psb_z_mat_mod class(psb_lz_base_sparse_mat), intent(inout) :: b end subroutine psb_lz_cp_to end interface + ! + ! Mixed type conversions + ! + interface + subroutine psb_lz_mv_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat + class(psb_lzspmat_type), intent(inout) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + end subroutine psb_lz_mv_from_ib + end interface + + interface + subroutine psb_lz_cp_from_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat + class(psb_lzspmat_type), intent(out) :: a + class(psb_z_base_sparse_mat), intent(in) :: b + end subroutine psb_lz_cp_from_ib + end interface + + interface + subroutine psb_lz_mv_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat + class(psb_lzspmat_type), intent(inout) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + end subroutine psb_lz_mv_to_ib + end interface + + interface + subroutine psb_lz_cp_to_ib(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat + class(psb_lzspmat_type), intent(in) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + end subroutine psb_lz_cp_to_ib + end interface + + interface + subroutine psb_lz_mv_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type + class(psb_lzspmat_type), intent(inout) :: a + class(psb_zspmat_type), intent(inout) :: b + end subroutine psb_lz_mv_from_i + end interface + + interface + subroutine psb_lz_cp_from_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type + class(psb_lzspmat_type), intent(out) :: a + class(psb_zspmat_type), intent(in) :: b + end subroutine psb_lz_cp_from_i + end interface + + interface + subroutine psb_lz_mv_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type + class(psb_lzspmat_type), intent(inout) :: a + class(psb_zspmat_type), intent(inout) :: b + end subroutine psb_lz_mv_to_i + end interface + + interface + subroutine psb_lz_cp_to_i(a,b) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type + class(psb_lzspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + end subroutine psb_lz_cp_to_i + end interface + ! ! Transfer the internal allocation to the target. diff --git a/base/modules/serial/psb_z_serial_mod.f90 b/base/modules/serial/psb_z_serial_mod.f90 index fb6d3dad9..c79740ab6 100644 --- a/base/modules/serial/psb_z_serial_mod.f90 +++ b/base/modules/serial/psb_z_serial_mod.f90 @@ -204,14 +204,14 @@ module psb_z_serial_mod end interface psb_aspxpby interface psb_spspmm -!!$ subroutine psb_zspspmm(a,b,c,info) -!!$ use psb_z_mat_mod, only : psb_zspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_zspmat_type), intent(in) :: a,b -!!$ type(psb_zspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_zspspmm + subroutine psb_lzspspmm(a,b,c,info) + use psb_z_mat_mod, only : psb_lzspmat_type + import :: psb_ipk_ + implicit none + type(psb_lzspmat_type), intent(in) :: a,b + type(psb_lzspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lzspspmm subroutine psb_lzcsrspspmm(a,b,c,info) use psb_z_mat_mod, only : psb_lz_csr_sparse_mat import :: psb_ipk_ @@ -231,14 +231,14 @@ module psb_z_serial_mod end interface psb_spspmm interface psb_symbmm -!!$ subroutine psb_zsymbmm(a,b,c,info) -!!$ use psb_z_mat_mod, only : psb_zspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_zspmat_type), intent(in) :: a,b -!!$ type(psb_zspmat_type), intent(out) :: c -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine psb_zsymbmm + subroutine psb_lzsymbmm(a,b,c,info) + use psb_z_mat_mod, only : psb_lzspmat_type + import :: psb_ipk_ + implicit none + type(psb_lzspmat_type), intent(in) :: a,b + type(psb_lzspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + end subroutine psb_lzsymbmm subroutine psb_lzbase_symbmm(a,b,c,info) use psb_z_mat_mod, only : psb_lz_base_sparse_mat, psb_lz_csr_sparse_mat import :: psb_ipk_ @@ -250,13 +250,13 @@ module psb_z_serial_mod end interface psb_symbmm interface psb_numbmm -!!$ subroutine psb_znumbmm(a,b,c) -!!$ use psb_z_mat_mod, only : psb_zspmat_type -!!$ import :: psb_ipk_ -!!$ implicit none -!!$ type(psb_zspmat_type), intent(in) :: a,b -!!$ type(psb_zspmat_type), intent(inout) :: c -!!$ end subroutine psb_znumbmm + subroutine psb_lznumbmm(a,b,c) + use psb_z_mat_mod, only : psb_lzspmat_type + import :: psb_ipk_ + implicit none + type(psb_lzspmat_type), intent(in) :: a,b + type(psb_lzspmat_type), intent(inout) :: c + end subroutine psb_lznumbmm subroutine psb_lzbase_numbmm(a,b,c) use psb_z_mat_mod, only : psb_lz_base_sparse_mat, psb_lz_csr_sparse_mat import :: psb_ipk_ diff --git a/base/serial/Makefile b/base/serial/Makefile index 7f3b7c990..1ce9156bb 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -8,7 +8,7 @@ FOBJS = psb_lsame.o psi_m_serial_impl.o psi_e_serial_impl.o \ psb_sspspmm.o psb_dspspmm.o psb_cspspmm.o psb_zspspmm.o \ psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \ psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \ - smmp.o \ + smmp.o lsmmp.o \ psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o\ psb_spdot_srtd.o psb_aspxpby.o psb_spge_dot.o\ psb_sgelp.o psb_dgelp.o psb_cgelp.o psb_zgelp.o \ diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index 7a1b18e70..b9a8d363d 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -2432,6 +2432,155 @@ subroutine psb_c_scals(d,a,info) end subroutine psb_c_scals +subroutine psb_c_mv_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_mv_from_lb + implicit none + + class(psb_cspmat_type), intent(inout) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_lfmt(b,info) + +end subroutine psb_c_mv_from_lb + + +subroutine psb_c_cp_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_cp_from_lb + implicit none + + class(psb_cspmat_type), intent(inout) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_lfmt(b,info) + +end subroutine psb_c_cp_from_lb + +subroutine psb_c_mv_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_mv_to_lb + implicit none + + class(psb_cspmat_type), intent(inout) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_lfmt(b,info) + call a%free() + end if + +end subroutine psb_c_mv_to_lb + +subroutine psb_c_cp_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_cp_to_lb + implicit none + class(psb_cspmat_type), intent(in) :: a + class(psb_lc_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_lfmt(b,info) + end if + +end subroutine psb_c_cp_to_lb + +subroutine psb_c_mv_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_mv_from_l + implicit none + class(psb_cspmat_type), intent(inout) :: a + class(psb_lcspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_lfmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_c_mv_from_l + + +subroutine psb_c_cp_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_cp_from_l + implicit none + + class(psb_cspmat_type), intent(out) :: a + class(psb_lcspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_lfmt(b%a,info) + else + call a%free() + end if +end subroutine psb_c_cp_from_l + +subroutine psb_c_mv_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_mv_to_l + implicit none + + class(psb_cspmat_type), intent(inout) :: a + class(psb_lcspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_lc_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_lfmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_c_mv_to_l + +subroutine psb_c_cp_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_c_cp_to_l + implicit none + + class(psb_cspmat_type), intent(in) :: a + class(psb_lcspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_lc_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_lfmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_c_cp_to_l + + +! +! +! lc versions +! + subroutine psb_lc_set_nrows(m,a) use psb_c_mat_mod, psb_protect_name => psb_lc_set_nrows @@ -4331,5 +4480,148 @@ subroutine psb_lc_scals(d,a,info) end subroutine psb_lc_scals +subroutine psb_lc_mv_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_mv_from_ib + implicit none + + class(psb_lcspmat_type), intent(inout) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_ifmt(b,info) + +end subroutine psb_lc_mv_from_ib + +subroutine psb_lc_cp_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_cp_from_ib + implicit none + + class(psb_lcspmat_type), intent(inout) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_ifmt(b,info) + +end subroutine psb_lc_cp_from_ib + +subroutine psb_lc_mv_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_mv_to_ib + implicit none + + class(psb_lcspmat_type), intent(inout) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_ifmt(b,info) + call a%free() + end if + +end subroutine psb_lc_mv_to_ib + +subroutine psb_lc_cp_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_cp_to_ib + implicit none + class(psb_lcspmat_type), intent(in) :: a + class(psb_c_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_ifmt(b,info) + end if + +end subroutine psb_lc_cp_to_ib + +subroutine psb_lc_mv_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_mv_from_i + implicit none + class(psb_lcspmat_type), intent(inout) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_ifmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_lc_mv_from_i + + +subroutine psb_lc_cp_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_cp_from_i + implicit none + + class(psb_lcspmat_type), intent(out) :: a + class(psb_cspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_ifmt(b%a,info) + else + call a%free() + end if +end subroutine psb_lc_cp_from_i + +subroutine psb_lc_mv_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_mv_to_i + implicit none + + class(psb_lcspmat_type), intent(inout) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_c_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_ifmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_lc_mv_to_i + +subroutine psb_lc_cp_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_c_mat_mod, psb_protect_name => psb_lc_cp_to_i + implicit none + + class(psb_lcspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_c_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_ifmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_lc_cp_to_i + + diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index a72dbb775..b8f54b1aa 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -2432,6 +2432,155 @@ subroutine psb_d_scals(d,a,info) end subroutine psb_d_scals +subroutine psb_d_mv_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_mv_from_lb + implicit none + + class(psb_dspmat_type), intent(inout) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_lfmt(b,info) + +end subroutine psb_d_mv_from_lb + + +subroutine psb_d_cp_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_cp_from_lb + implicit none + + class(psb_dspmat_type), intent(inout) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_lfmt(b,info) + +end subroutine psb_d_cp_from_lb + +subroutine psb_d_mv_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_mv_to_lb + implicit none + + class(psb_dspmat_type), intent(inout) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_lfmt(b,info) + call a%free() + end if + +end subroutine psb_d_mv_to_lb + +subroutine psb_d_cp_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_cp_to_lb + implicit none + class(psb_dspmat_type), intent(in) :: a + class(psb_ld_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_lfmt(b,info) + end if + +end subroutine psb_d_cp_to_lb + +subroutine psb_d_mv_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_mv_from_l + implicit none + class(psb_dspmat_type), intent(inout) :: a + class(psb_ldspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_lfmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_d_mv_from_l + + +subroutine psb_d_cp_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_cp_from_l + implicit none + + class(psb_dspmat_type), intent(out) :: a + class(psb_ldspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_lfmt(b%a,info) + else + call a%free() + end if +end subroutine psb_d_cp_from_l + +subroutine psb_d_mv_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_mv_to_l + implicit none + + class(psb_dspmat_type), intent(inout) :: a + class(psb_ldspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_ld_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_lfmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_d_mv_to_l + +subroutine psb_d_cp_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_d_cp_to_l + implicit none + + class(psb_dspmat_type), intent(in) :: a + class(psb_ldspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_ld_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_lfmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_d_cp_to_l + + +! +! +! ld versions +! + subroutine psb_ld_set_nrows(m,a) use psb_d_mat_mod, psb_protect_name => psb_ld_set_nrows @@ -4331,5 +4480,148 @@ subroutine psb_ld_scals(d,a,info) end subroutine psb_ld_scals +subroutine psb_ld_mv_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_mv_from_ib + implicit none + + class(psb_ldspmat_type), intent(inout) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_ifmt(b,info) + +end subroutine psb_ld_mv_from_ib + +subroutine psb_ld_cp_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_cp_from_ib + implicit none + + class(psb_ldspmat_type), intent(inout) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_ifmt(b,info) + +end subroutine psb_ld_cp_from_ib + +subroutine psb_ld_mv_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_mv_to_ib + implicit none + + class(psb_ldspmat_type), intent(inout) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_ifmt(b,info) + call a%free() + end if + +end subroutine psb_ld_mv_to_ib + +subroutine psb_ld_cp_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_cp_to_ib + implicit none + class(psb_ldspmat_type), intent(in) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_ifmt(b,info) + end if + +end subroutine psb_ld_cp_to_ib + +subroutine psb_ld_mv_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_mv_from_i + implicit none + class(psb_ldspmat_type), intent(inout) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_ifmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_ld_mv_from_i + + +subroutine psb_ld_cp_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_cp_from_i + implicit none + + class(psb_ldspmat_type), intent(out) :: a + class(psb_dspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_ifmt(b%a,info) + else + call a%free() + end if +end subroutine psb_ld_cp_from_i + +subroutine psb_ld_mv_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_mv_to_i + implicit none + + class(psb_ldspmat_type), intent(inout) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_d_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_ifmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_ld_mv_to_i + +subroutine psb_ld_cp_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_d_mat_mod, psb_protect_name => psb_ld_cp_to_i + implicit none + + class(psb_ldspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_d_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_ifmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_ld_cp_to_i + + diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index 5907f9c0b..af9235641 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -2432,6 +2432,155 @@ subroutine psb_s_scals(d,a,info) end subroutine psb_s_scals +subroutine psb_s_mv_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_mv_from_lb + implicit none + + class(psb_sspmat_type), intent(inout) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_lfmt(b,info) + +end subroutine psb_s_mv_from_lb + + +subroutine psb_s_cp_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_cp_from_lb + implicit none + + class(psb_sspmat_type), intent(inout) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_lfmt(b,info) + +end subroutine psb_s_cp_from_lb + +subroutine psb_s_mv_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_mv_to_lb + implicit none + + class(psb_sspmat_type), intent(inout) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_lfmt(b,info) + call a%free() + end if + +end subroutine psb_s_mv_to_lb + +subroutine psb_s_cp_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_cp_to_lb + implicit none + class(psb_sspmat_type), intent(in) :: a + class(psb_ls_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_lfmt(b,info) + end if + +end subroutine psb_s_cp_to_lb + +subroutine psb_s_mv_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_mv_from_l + implicit none + class(psb_sspmat_type), intent(inout) :: a + class(psb_lsspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_lfmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_s_mv_from_l + + +subroutine psb_s_cp_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_cp_from_l + implicit none + + class(psb_sspmat_type), intent(out) :: a + class(psb_lsspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_lfmt(b%a,info) + else + call a%free() + end if +end subroutine psb_s_cp_from_l + +subroutine psb_s_mv_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_mv_to_l + implicit none + + class(psb_sspmat_type), intent(inout) :: a + class(psb_lsspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_ls_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_lfmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_s_mv_to_l + +subroutine psb_s_cp_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_s_cp_to_l + implicit none + + class(psb_sspmat_type), intent(in) :: a + class(psb_lsspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_ls_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_lfmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_s_cp_to_l + + +! +! +! ls versions +! + subroutine psb_ls_set_nrows(m,a) use psb_s_mat_mod, psb_protect_name => psb_ls_set_nrows @@ -4331,5 +4480,148 @@ subroutine psb_ls_scals(d,a,info) end subroutine psb_ls_scals +subroutine psb_ls_mv_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_mv_from_ib + implicit none + + class(psb_lsspmat_type), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_ifmt(b,info) + +end subroutine psb_ls_mv_from_ib + +subroutine psb_ls_cp_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_cp_from_ib + implicit none + + class(psb_lsspmat_type), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_ifmt(b,info) + +end subroutine psb_ls_cp_from_ib + +subroutine psb_ls_mv_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_mv_to_ib + implicit none + + class(psb_lsspmat_type), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_ifmt(b,info) + call a%free() + end if + +end subroutine psb_ls_mv_to_ib + +subroutine psb_ls_cp_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_cp_to_ib + implicit none + class(psb_lsspmat_type), intent(in) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_ifmt(b,info) + end if + +end subroutine psb_ls_cp_to_ib + +subroutine psb_ls_mv_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_mv_from_i + implicit none + class(psb_lsspmat_type), intent(inout) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_ifmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_ls_mv_from_i + + +subroutine psb_ls_cp_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_cp_from_i + implicit none + + class(psb_lsspmat_type), intent(out) :: a + class(psb_sspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_ifmt(b%a,info) + else + call a%free() + end if +end subroutine psb_ls_cp_from_i + +subroutine psb_ls_mv_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_mv_to_i + implicit none + + class(psb_lsspmat_type), intent(inout) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_s_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_ifmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_ls_mv_to_i + +subroutine psb_ls_cp_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_s_mat_mod, psb_protect_name => psb_ls_cp_to_i + implicit none + + class(psb_lsspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_s_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_ifmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_ls_cp_to_i + + diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index 37e6db48e..bf4b266db 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -2432,6 +2432,155 @@ subroutine psb_z_scals(d,a,info) end subroutine psb_z_scals +subroutine psb_z_mv_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_mv_from_lb + implicit none + + class(psb_zspmat_type), intent(inout) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_lfmt(b,info) + +end subroutine psb_z_mv_from_lb + + +subroutine psb_z_cp_from_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_cp_from_lb + implicit none + + class(psb_zspmat_type), intent(inout) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_lfmt(b,info) + +end subroutine psb_z_cp_from_lb + +subroutine psb_z_mv_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_mv_to_lb + implicit none + + class(psb_zspmat_type), intent(inout) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_lfmt(b,info) + call a%free() + end if + +end subroutine psb_z_mv_to_lb + +subroutine psb_z_cp_to_lb(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_cp_to_lb + implicit none + class(psb_zspmat_type), intent(in) :: a + class(psb_lz_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_lfmt(b,info) + end if + +end subroutine psb_z_cp_to_lb + +subroutine psb_z_mv_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_mv_from_l + implicit none + class(psb_zspmat_type), intent(inout) :: a + class(psb_lzspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_lfmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_z_mv_from_l + + +subroutine psb_z_cp_from_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_cp_from_l + implicit none + + class(psb_zspmat_type), intent(out) :: a + class(psb_lzspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_lfmt(b%a,info) + else + call a%free() + end if +end subroutine psb_z_cp_from_l + +subroutine psb_z_mv_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_mv_to_l + implicit none + + class(psb_zspmat_type), intent(inout) :: a + class(psb_lzspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_lz_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_lfmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_z_mv_to_l + +subroutine psb_z_cp_to_l(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_z_cp_to_l + implicit none + + class(psb_zspmat_type), intent(in) :: a + class(psb_lzspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_lz_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_lfmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_z_cp_to_l + + +! +! +! lz versions +! + subroutine psb_lz_set_nrows(m,a) use psb_z_mat_mod, psb_protect_name => psb_lz_set_nrows @@ -4331,5 +4480,148 @@ subroutine psb_lz_scals(d,a,info) end subroutine psb_lz_scals +subroutine psb_lz_mv_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_mv_from_ib + implicit none + + class(psb_lzspmat_type), intent(inout) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%mv_from_ifmt(b,info) + +end subroutine psb_lz_mv_from_ib + +subroutine psb_lz_cp_from_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_cp_from_ib + implicit none + + class(psb_lzspmat_type), intent(inout) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info) + if (info == psb_success_) call a%a%cp_from_ifmt(b,info) + +end subroutine psb_lz_cp_from_ib + +subroutine psb_lz_mv_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_mv_to_ib + implicit none + + class(psb_lzspmat_type), intent(inout) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%mv_to_ifmt(b,info) + call a%free() + end if + +end subroutine psb_lz_mv_to_ib + +subroutine psb_lz_cp_to_ib(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_cp_to_ib + implicit none + class(psb_lzspmat_type), intent(in) :: a + class(psb_z_base_sparse_mat), intent(inout) :: b + integer(psb_ipk_) :: info + + if (.not.allocated(a%a)) then + call b%free() + else + call a%a%cp_to_ifmt(b,info) + end if + +end subroutine psb_lz_cp_to_ib + +subroutine psb_lz_mv_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_mv_from_i + implicit none + class(psb_lzspmat_type), intent(inout) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info) + call a%a%mv_from_ifmt(b%a,info) + else + call a%free() + end if + call b%free() + +end subroutine psb_lz_mv_from_i + + +subroutine psb_lz_cp_from_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_cp_from_i + implicit none + + class(psb_lzspmat_type), intent(out) :: a + class(psb_zspmat_type), intent(in) :: b + integer(psb_ipk_) :: info + + if (allocated(b%a)) then + if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info) + call a%a%cp_from_ifmt(b%a,info) + else + call a%free() + end if +end subroutine psb_lz_cp_from_i + +subroutine psb_lz_mv_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_mv_to_i + implicit none + + class(psb_lzspmat_type), intent(inout) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_z_csr_sparse_mat :: b%a, stat=info) + call a%a%mv_to_ifmt(b%a,info) + else + call b%free() + end if + call a%free() + +end subroutine psb_lz_mv_to_i + +subroutine psb_lz_cp_to_i(a,b) + use psb_error_mod + use psb_const_mod + use psb_z_mat_mod, psb_protect_name => psb_lz_cp_to_i + implicit none + + class(psb_lzspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_) :: info + + if (allocated(a%a)) then + if (.not.allocated(b%a)) allocate(psb_z_csr_sparse_mat :: b%a, stat=info) + call a%a%cp_to_ifmt(b%a,info) + else + call b%free() + end if + +end subroutine psb_lz_cp_to_i + + diff --git a/base/serial/psb_cnumbmm.f90 b/base/serial/psb_cnumbmm.f90 index 94e6e1a5a..c965d4f3a 100644 --- a/base/serial/psb_cnumbmm.f90 +++ b/base/serial/psb_cnumbmm.f90 @@ -39,7 +39,6 @@ ! rewritten in Fortran 95/2003 making use of our sparse matrix facilities. ! ! - subroutine psb_cnumbmm(a,b,c) use psb_base_mod, psb_protect_name => psb_cnumbmm implicit none @@ -234,3 +233,200 @@ contains end subroutine gen_numbmm end subroutine psb_cbase_numbmm + + + +subroutine psb_lcnumbmm(a,b,c) + use psb_base_mod, psb_protect_name => psb_lcnumbmm + implicit none + + type(psb_lcspmat_type), intent(in) :: a,b + type(psb_lcspmat_type), intent(inout) :: c + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_lc_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lcnumbmm + +subroutine psb_lcbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_lcbase_numbmm + implicit none + + class(psb_lc_base_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + complex(psb_spk_), allocatable :: temp(:) + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + endif + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + + ! + ! Note: we still have to test about possible performance hits. + ! + ! + call psb_ensure_size(ione*size(c%ja),c%val,info) + select type(a) + type is (psb_lc_csr_sparse_mat) + select type(b) + type is (psb_lc_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call c%set_asb() + deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_lc_csr_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(inout) :: c + complex(psb_spk_) :: temp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call lcnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,& + & b%irp,b%ja,lzero,b%val,& + & c%irp,c%ja,lzero,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_lc_base_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_) :: info + complex(psb_spk_) :: temp(:) + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + complex(psb_spk_), allocatable :: aval(:),bval(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln + complex(psb_spk_) :: ajj + + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() + maxlmn = max(l,m,n) + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & aval(maxlmn),bval(maxlmn), stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i = 1,maxlmn + temp(i) = czero + end do + minlm = min(l,m) + minln = min(l,n) + minmn = min(m,n) + do i = 1,n + + call a%csget(i,i,nazr,iarw,iacl,aval,info) + do jj=1, nazr + j=iacl(jj) + ajj = aval(jj) + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + + endif + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = psb_err_pivot_too_small_ + return + else + temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) + endif + enddo + end do + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = psb_err_invalid_ovr_num_ + return + else + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = czero + endif + end do + end do + + + end subroutine gen_numbmm + +end subroutine psb_lcbase_numbmm diff --git a/base/serial/psb_csymbmm.f90 b/base/serial/psb_csymbmm.f90 index 3343c61de..25818791d 100644 --- a/base/serial/psb_csymbmm.f90 +++ b/base/serial/psb_csymbmm.f90 @@ -255,3 +255,223 @@ contains end subroutine gen_symbmm end subroutine psb_cbase_symbmm + + + +subroutine psb_lcsymbmm(a,b,c,info) + use psb_base_mod, psb_protect_name => psb_lcsymbmm + implicit none + + type(psb_lcspmat_type), intent(in) :: a,b + type(psb_lcspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + type(psb_lc_csr_sparse_mat), allocatable :: ccsr + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr,stat=info) + + if (info == psb_success_) then + call psb_symbmm(a%a,b%a,ccsr,info) + else + info = psb_err_alloc_dealloc_ + end if + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lcsymbmm + +subroutine psb_lcbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_lcbase_symbmm + implicit none + + class(psb_lc_base_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + integer(psb_ipk_) :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + info = psb_err_invalid_matrix_sizes_ + call psb_errpush(info,name) + goto 9999 + endif + allocate(itemp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + ! + ! Note: we need to test whether there is a performance impact + ! in not using the original Douglas & Bank code. + ! + select type(a) + type is (psb_lc_csr_sparse_mat) + select type(b) + type is (psb_lc_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) + deallocate(itemp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_lc_csr_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: itemp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,& + & b%irp,b%ja,lzero,& + & c%irp,c%ja,lzero,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_lc_base_sparse_mat), intent(in) :: a,b + type(psb_lc_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: index(:) + integer(psb_ipk_) :: info + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer(psb_lpk_) :: nze, ma,na,mb,nb + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb + maxlmn = max(l,m,n) + + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i=1,maxlmn + index(i)=0 + end do + + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) + + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=psb_err_pivot_too_small_ + return + else + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif + end do + end do + + c%irp(i+1)=c%irp(i)+length + + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_lcbase_symbmm diff --git a/base/serial/psb_dnumbmm.f90 b/base/serial/psb_dnumbmm.f90 index 2de4b9039..c1d3951cd 100644 --- a/base/serial/psb_dnumbmm.f90 +++ b/base/serial/psb_dnumbmm.f90 @@ -39,7 +39,6 @@ ! rewritten in Fortran 95/2003 making use of our sparse matrix facilities. ! ! - subroutine psb_dnumbmm(a,b,c) use psb_base_mod, psb_protect_name => psb_dnumbmm implicit none @@ -234,3 +233,200 @@ contains end subroutine gen_numbmm end subroutine psb_dbase_numbmm + + + +subroutine psb_ldnumbmm(a,b,c) + use psb_base_mod, psb_protect_name => psb_ldnumbmm + implicit none + + type(psb_ldspmat_type), intent(in) :: a,b + type(psb_ldspmat_type), intent(inout) :: c + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_ld_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_ldnumbmm + +subroutine psb_ldbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_ldbase_numbmm + implicit none + + class(psb_ld_base_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + real(psb_dpk_), allocatable :: temp(:) + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + endif + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + + ! + ! Note: we still have to test about possible performance hits. + ! + ! + call psb_ensure_size(ione*size(c%ja),c%val,info) + select type(a) + type is (psb_ld_csr_sparse_mat) + select type(b) + type is (psb_ld_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call c%set_asb() + deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_ld_csr_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(inout) :: c + real(psb_dpk_) :: temp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call ldnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,& + & b%irp,b%ja,lzero,b%val,& + & c%irp,c%ja,lzero,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_ld_base_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_) :: info + real(psb_dpk_) :: temp(:) + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + real(psb_dpk_), allocatable :: aval(:),bval(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln + real(psb_dpk_) :: ajj + + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() + maxlmn = max(l,m,n) + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & aval(maxlmn),bval(maxlmn), stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i = 1,maxlmn + temp(i) = dzero + end do + minlm = min(l,m) + minln = min(l,n) + minmn = min(m,n) + do i = 1,n + + call a%csget(i,i,nazr,iarw,iacl,aval,info) + do jj=1, nazr + j=iacl(jj) + ajj = aval(jj) + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + + endif + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = psb_err_pivot_too_small_ + return + else + temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) + endif + enddo + end do + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = psb_err_invalid_ovr_num_ + return + else + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = dzero + endif + end do + end do + + + end subroutine gen_numbmm + +end subroutine psb_ldbase_numbmm diff --git a/base/serial/psb_dsymbmm.f90 b/base/serial/psb_dsymbmm.f90 index 848a5cfd0..3dcad00f9 100644 --- a/base/serial/psb_dsymbmm.f90 +++ b/base/serial/psb_dsymbmm.f90 @@ -255,3 +255,223 @@ contains end subroutine gen_symbmm end subroutine psb_dbase_symbmm + + + +subroutine psb_ldsymbmm(a,b,c,info) + use psb_base_mod, psb_protect_name => psb_ldsymbmm + implicit none + + type(psb_ldspmat_type), intent(in) :: a,b + type(psb_ldspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + type(psb_ld_csr_sparse_mat), allocatable :: ccsr + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr,stat=info) + + if (info == psb_success_) then + call psb_symbmm(a%a,b%a,ccsr,info) + else + info = psb_err_alloc_dealloc_ + end if + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_ldsymbmm + +subroutine psb_ldbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_ldbase_symbmm + implicit none + + class(psb_ld_base_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + integer(psb_ipk_) :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + info = psb_err_invalid_matrix_sizes_ + call psb_errpush(info,name) + goto 9999 + endif + allocate(itemp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + ! + ! Note: we need to test whether there is a performance impact + ! in not using the original Douglas & Bank code. + ! + select type(a) + type is (psb_ld_csr_sparse_mat) + select type(b) + type is (psb_ld_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) + deallocate(itemp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_ld_csr_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: itemp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,& + & b%irp,b%ja,lzero,& + & c%irp,c%ja,lzero,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_ld_base_sparse_mat), intent(in) :: a,b + type(psb_ld_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: index(:) + integer(psb_ipk_) :: info + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer(psb_lpk_) :: nze, ma,na,mb,nb + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb + maxlmn = max(l,m,n) + + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i=1,maxlmn + index(i)=0 + end do + + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) + + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=psb_err_pivot_too_small_ + return + else + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif + end do + end do + + c%irp(i+1)=c%irp(i)+length + + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_ldbase_symbmm diff --git a/base/serial/psb_snumbmm.f90 b/base/serial/psb_snumbmm.f90 index 9b37a10af..ceffb977c 100644 --- a/base/serial/psb_snumbmm.f90 +++ b/base/serial/psb_snumbmm.f90 @@ -39,7 +39,6 @@ ! rewritten in Fortran 95/2003 making use of our sparse matrix facilities. ! ! - subroutine psb_snumbmm(a,b,c) use psb_base_mod, psb_protect_name => psb_snumbmm implicit none @@ -234,3 +233,200 @@ contains end subroutine gen_numbmm end subroutine psb_sbase_numbmm + + + +subroutine psb_lsnumbmm(a,b,c) + use psb_base_mod, psb_protect_name => psb_lsnumbmm + implicit none + + type(psb_lsspmat_type), intent(in) :: a,b + type(psb_lsspmat_type), intent(inout) :: c + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_ls_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lsnumbmm + +subroutine psb_lsbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_lsbase_numbmm + implicit none + + class(psb_ls_base_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + real(psb_spk_), allocatable :: temp(:) + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + endif + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + + ! + ! Note: we still have to test about possible performance hits. + ! + ! + call psb_ensure_size(ione*size(c%ja),c%val,info) + select type(a) + type is (psb_ls_csr_sparse_mat) + select type(b) + type is (psb_ls_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call c%set_asb() + deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_ls_csr_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(inout) :: c + real(psb_spk_) :: temp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call lsnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,& + & b%irp,b%ja,lzero,b%val,& + & c%irp,c%ja,lzero,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_ls_base_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_) :: info + real(psb_spk_) :: temp(:) + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + real(psb_spk_), allocatable :: aval(:),bval(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln + real(psb_spk_) :: ajj + + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() + maxlmn = max(l,m,n) + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & aval(maxlmn),bval(maxlmn), stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i = 1,maxlmn + temp(i) = szero + end do + minlm = min(l,m) + minln = min(l,n) + minmn = min(m,n) + do i = 1,n + + call a%csget(i,i,nazr,iarw,iacl,aval,info) + do jj=1, nazr + j=iacl(jj) + ajj = aval(jj) + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + + endif + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = psb_err_pivot_too_small_ + return + else + temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) + endif + enddo + end do + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = psb_err_invalid_ovr_num_ + return + else + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = szero + endif + end do + end do + + + end subroutine gen_numbmm + +end subroutine psb_lsbase_numbmm diff --git a/base/serial/psb_ssymbmm.f90 b/base/serial/psb_ssymbmm.f90 index e9d10c0b1..729dd856b 100644 --- a/base/serial/psb_ssymbmm.f90 +++ b/base/serial/psb_ssymbmm.f90 @@ -255,3 +255,223 @@ contains end subroutine gen_symbmm end subroutine psb_sbase_symbmm + + + +subroutine psb_lssymbmm(a,b,c,info) + use psb_base_mod, psb_protect_name => psb_lssymbmm + implicit none + + type(psb_lsspmat_type), intent(in) :: a,b + type(psb_lsspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + type(psb_ls_csr_sparse_mat), allocatable :: ccsr + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr,stat=info) + + if (info == psb_success_) then + call psb_symbmm(a%a,b%a,ccsr,info) + else + info = psb_err_alloc_dealloc_ + end if + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lssymbmm + +subroutine psb_lsbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_lsbase_symbmm + implicit none + + class(psb_ls_base_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + integer(psb_ipk_) :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + info = psb_err_invalid_matrix_sizes_ + call psb_errpush(info,name) + goto 9999 + endif + allocate(itemp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + ! + ! Note: we need to test whether there is a performance impact + ! in not using the original Douglas & Bank code. + ! + select type(a) + type is (psb_ls_csr_sparse_mat) + select type(b) + type is (psb_ls_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) + deallocate(itemp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_ls_csr_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: itemp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,& + & b%irp,b%ja,lzero,& + & c%irp,c%ja,lzero,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_ls_base_sparse_mat), intent(in) :: a,b + type(psb_ls_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: index(:) + integer(psb_ipk_) :: info + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer(psb_lpk_) :: nze, ma,na,mb,nb + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb + maxlmn = max(l,m,n) + + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i=1,maxlmn + index(i)=0 + end do + + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) + + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=psb_err_pivot_too_small_ + return + else + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif + end do + end do + + c%irp(i+1)=c%irp(i)+length + + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_lsbase_symbmm diff --git a/base/serial/psb_znumbmm.f90 b/base/serial/psb_znumbmm.f90 index 130de6ebe..be4e10263 100644 --- a/base/serial/psb_znumbmm.f90 +++ b/base/serial/psb_znumbmm.f90 @@ -39,7 +39,6 @@ ! rewritten in Fortran 95/2003 making use of our sparse matrix facilities. ! ! - subroutine psb_znumbmm(a,b,c) use psb_base_mod, psb_protect_name => psb_znumbmm implicit none @@ -234,3 +233,200 @@ contains end subroutine gen_numbmm end subroutine psb_zbase_numbmm + + + +subroutine psb_lznumbmm(a,b,c) + use psb_base_mod, psb_protect_name => psb_lznumbmm + implicit none + + type(psb_lzspmat_type), intent(in) :: a,b + type(psb_lzspmat_type), intent(inout) :: c + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_numbmm' + + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + select type(aa=>c%a) + type is (psb_lz_csr_sparse_mat) + call psb_numbmm(a%a,b%a,aa) + class default + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end select + + call c%set_asb() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lznumbmm + +subroutine psb_lzbase_numbmm(a,b,c) + use psb_mat_mod + use psb_string_mod + use psb_serial_mod, psb_protect_name => psb_lzbase_numbmm + implicit none + + class(psb_lz_base_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + complex(psb_dpk_), allocatable :: temp(:) + integer(psb_ipk_) :: info + integer(psb_ipk_) :: err_act + name='psb_numbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + endif + allocate(temp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + + ! + ! Note: we still have to test about possible performance hits. + ! + ! + call psb_ensure_size(ione*size(c%ja),c%val,info) + select type(a) + type is (psb_lz_csr_sparse_mat) + select type(b) + type is (psb_lz_csr_sparse_mat) + call csr_numbmm(a,b,c,temp,info) + class default + call gen_numbmm(a,b,c,temp,info) + end select + class default + call gen_numbmm(a,b,c,temp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call c%set_asb() + deallocate(temp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_numbmm(a,b,c,temp,info) + type(psb_lz_csr_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(inout) :: c + complex(psb_dpk_) :: temp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + call lznumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,& + & b%irp,b%ja,lzero,b%val,& + & c%irp,c%ja,lzero,c%val,temp) + + + end subroutine csr_numbmm + + subroutine gen_numbmm(a,b,c,temp,info) + class(psb_lz_base_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(inout) :: c + integer(psb_ipk_) :: info + complex(psb_dpk_) :: temp(:) + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + complex(psb_dpk_), allocatable :: aval(:),bval(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln + complex(psb_dpk_) :: ajj + + n = a%get_nrows() + m = a%get_ncols() + l = b%get_ncols() + maxlmn = max(l,m,n) + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & aval(maxlmn),bval(maxlmn), stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i = 1,maxlmn + temp(i) = zzero + end do + minlm = min(l,m) + minln = min(l,n) + minmn = min(m,n) + do i = 1,n + + call a%csget(i,i,nazr,iarw,iacl,aval,info) + do jj=1, nazr + j=iacl(jj) + ajj = aval(jj) + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m + info = 1 + return + + endif + call b%csget(j,j,nbzr,ibrw,ibcl,bval,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn + info = psb_err_pivot_too_small_ + return + else + temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k) + endif + enddo + end do + do j = c%irp(i),c%irp(i+1)-1 + if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then + write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn + info = psb_err_invalid_ovr_num_ + return + else + c%val(j) = temp(c%ja(j)) + temp(c%ja(j)) = zzero + endif + end do + end do + + + end subroutine gen_numbmm + +end subroutine psb_lzbase_numbmm diff --git a/base/serial/psb_zsymbmm.f90 b/base/serial/psb_zsymbmm.f90 index 67094aaf4..ada823264 100644 --- a/base/serial/psb_zsymbmm.f90 +++ b/base/serial/psb_zsymbmm.f90 @@ -255,3 +255,223 @@ contains end subroutine gen_symbmm end subroutine psb_zbase_symbmm + + + +subroutine psb_lzsymbmm(a,b,c,info) + use psb_base_mod, psb_protect_name => psb_lzsymbmm + implicit none + + type(psb_lzspmat_type), intent(in) :: a,b + type(psb_lzspmat_type), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + type(psb_lz_csr_sparse_mat), allocatable :: ccsr + integer(psb_ipk_) :: err_act + character(len=*), parameter :: name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + if ((a%is_null()) .or.(b%is_null())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(ccsr,stat=info) + + if (info == psb_success_) then + call psb_symbmm(a%a,b%a,ccsr,info) + else + info = psb_err_alloc_dealloc_ + end if + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + call move_alloc(ccsr,c%a) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_lzsymbmm + +subroutine psb_lzbase_symbmm(a,b,c,info) + use psb_mat_mod + use psb_serial_mod, psb_protect_name => psb_lzbase_symbmm + implicit none + + class(psb_lz_base_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(out) :: c + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), allocatable :: itemp(:) + integer(psb_lpk_) :: nze, ma,na,mb,nb + character(len=20) :: name + integer(psb_ipk_) :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + info = psb_success_ + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + + if ( mb /= na ) then + write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb + info = psb_err_invalid_matrix_sizes_ + call psb_errpush(info,name) + goto 9999 + endif + allocate(itemp(max(ma,na,mb,nb)),stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_Errpush(info,name) + goto 9999 + endif + ! + ! Note: we need to test whether there is a performance impact + ! in not using the original Douglas & Bank code. + ! + select type(a) + type is (psb_lz_csr_sparse_mat) + select type(b) + type is (psb_lz_csr_sparse_mat) + call csr_symbmm(a,b,c,itemp,info) + class default + call gen_symbmm(a,b,c,itemp,info) + end select + class default + call gen_symbmm(a,b,c,itemp,info) + end select + + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_realloc(size(c%ja),c%val,info) + deallocate(itemp) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +contains + + subroutine csr_symbmm(a,b,c,itemp,info) + type(psb_lz_csr_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: itemp(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_) :: nze, ma,na,mb,nb + + info = psb_success_ + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,& + & b%irp,b%ja,lzero,& + & c%irp,c%ja,lzero,itemp) + + end subroutine csr_symbmm + subroutine gen_symbmm(a,b,c,index,info) + class(psb_lz_base_sparse_mat), intent(in) :: a,b + type(psb_lz_csr_sparse_mat), intent(out) :: c + integer(psb_lpk_) :: index(:) + integer(psb_ipk_) :: info + integer(psb_lpk_), allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:) + integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn + integer(psb_lpk_) :: nze, ma,na,mb,nb + + ma = a%get_nrows() + na = a%get_ncols() + mb = b%get_nrows() + nb = b%get_ncols() + + nze = max(ma+1,2*ma) + call c%allocate(ma,nb,nze) + + n = ma + m = na + l = nb + maxlmn = max(l,m,n) + + allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),& + & stat=info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + return + endif + + do i=1,maxlmn + index(i)=0 + end do + + c%irp(1)=1 + minlm = min(l,m) + minmn = min(m,n) + + main: do i=1,n + istart=-1 + length=0 + call a%csget(i,i,nazr,iarw,iacl,info) + do jj=1, nazr + + j=iacl(jj) + + if ((j<1).or.(j>m)) then + write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m + info = 1 + return + endif + call b%csget(j,j,nbzr,ibrw,ibcl,info) + do k=1,nbzr + if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then + write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn + info=psb_err_pivot_too_small_ + return + else + if(index(ibcl(k)) == 0) then + index(ibcl(k))=istart + istart=ibcl(k) + length=length+1 + endif + endif + end do + end do + + c%irp(i+1)=c%irp(i)+length + + if (c%irp(i+1) > size(c%ja)) then + if (n > (2*i)) then + nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i)) + else + nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) ) + endif + call psb_realloc(nze,c%ja,info) + end if + do j= c%irp(i),c%irp(i+1)-1 + c%ja(j)=istart + istart=index(istart) + index(c%ja(j))=0 + end do + call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1)) + index(i) = 0 + end do main + + end subroutine gen_symbmm + +end subroutine psb_lzbase_symbmm