diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 59788756..39109db1 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -71,8 +71,6 @@ module psbn_d_base_mat_mod procedure, pass(a) :: free => d_coo_free procedure, pass(a) :: trim => d_coo_trim procedure, pass(a) :: d_csgetrow => d_coo_csgetrow - procedure, pass(a) :: d_csgetblk => d_coo_csgetblk - procedure, pass(a) :: csclip => d_coo_csclip procedure, pass(a) :: print => d_coo_print procedure, pass(a) :: get_fmt => d_coo_get_fmt @@ -84,8 +82,7 @@ module psbn_d_base_mat_mod & d_fix_coo, d_coo_free, d_coo_print, d_coo_get_fmt, & & d_cp_coo_to_coo, d_cp_coo_from_coo, & & d_cp_coo_to_fmt, d_cp_coo_from_fmt, & - & d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_csgetblk, & - & d_coo_csclip + & d_coo_scals, d_coo_scal, d_coo_csgetrow interface @@ -568,6 +565,7 @@ contains end subroutine d_csgetrow + subroutine d_csgetblk(imin,imax,a,b,info,& & jmin,jmax,iren,append,rscale,cscale) ! Output is always in COO format @@ -583,19 +581,44 @@ contains integer, intent(in), optional :: iren(:) integer, intent(in), optional :: jmin,jmax logical, intent(in), optional :: rscale,cscale - Integer :: err_act + Integer :: err_act, nzin, nzout character(len=20) :: name='csget' + logical :: append_ logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 - call psb_get_erraction(err_act) - ! This is the base version. If we get here - ! it means the derived class is incomplete, - ! so we throw an error. - info = 700 - call psb_errpush(info,name,a_err=a%get_fmt()) + if (present(append)) then + append_ = append + else + append_ = .false. + endif + if (append_) then + nzin = a%get_nzeros() + else + nzin = 0 + endif - if (err_act /= psb_act_ret_) then + call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin, jmax=jmax, iren=iren, append=append_, & + & nzin=nzin, rscale=rscale, cscale=cscale) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then call psb_error() + return end if return @@ -614,19 +637,78 @@ contains integer,intent(out) :: info integer, intent(in), optional :: imin,imax,jmin,jmax logical, intent(in), optional :: rscale,cscale - Integer :: err_act - character(len=20) :: name='csclip' + + Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb + character(len=20) :: name='csget' + logical :: rscale_, cscale_ logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 - call psb_get_erraction(err_act) - ! This is the base version. If we get here - ! it means the derived class is incomplete, - ! so we throw an error. - info = 700 - call psb_errpush(info,name,a_err=a%get_fmt()) + nzin = 0 + if (present(imin)) then + imin_ = imin + else + imin_ = 1 + end if + if (present(imax)) then + imax_ = imax + else + imax_ = a%get_nrows() + end if + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + end if + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + end if + if (present(rscale)) then + rscale_ = rscale + else + rscale_ = .true. + end if + if (present(cscale)) then + cscale_ = cscale + else + cscale_ = .true. + end if - if (err_act /= psb_act_ret_) then + if (rscale_) then + mb = imax_ - imin_ +1 + else + mb = a%get_nrows() ! Should this be imax_ ?? + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = a%get_ncols() ! Should this be jmax_ ?? + endif + call b%allocate(mb,nb) + + call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=jmax_, append=.false., & + & nzin=nzin, rscale=rscale_, cscale=cscale_) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + + if (info /= 0) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then call psb_error() + return end if return @@ -1392,156 +1474,6 @@ contains end subroutine d_coo_csgetrow - subroutine d_coo_csgetblk(imin,imax,a,b,info,& - & jmin,jmax,iren,append,rscale,cscale) - ! Output is always in COO format - use psb_error_mod - use psb_const_mod - implicit none - - class(psbn_d_coo_sparse_mat), intent(in) :: a - class(psbn_d_coo_sparse_mat), intent(inout) :: b - integer, intent(in) :: imin,imax - integer,intent(out) :: info - logical, intent(in), optional :: append - integer, intent(in), optional :: iren(:) - integer, intent(in), optional :: jmin,jmax - logical, intent(in), optional :: rscale,cscale - Integer :: err_act, nzin, nzout - character(len=20) :: name='csget' - logical :: append_ - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - info = 0 - - if (present(append)) then - append_ = append - else - append_ = .false. - endif - if (append_) then - nzin = a%get_nzeros() - else - nzin = 0 - endif - - call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin, jmax=jmax, iren=iren, append=append_, & - & nzin=nzin, rscale=rscale, cscale=cscale) - - if (info /= 0) goto 9999 - - call b%set_nzeros(nzin+nzout) - call b%fix(info) - if (info /= 0) goto 9999 - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_coo_csgetblk - - - subroutine d_coo_csclip(a,b,info,& - & imin,imax,jmin,jmax,rscale,cscale) - ! Output is always in COO format - use psb_error_mod - use psb_const_mod - implicit none - - class(psbn_d_coo_sparse_mat), intent(in) :: a - class(psbn_d_coo_sparse_mat), intent(out) :: b - integer,intent(out) :: info - integer, intent(in), optional :: imin,imax,jmin,jmax - logical, intent(in), optional :: rscale,cscale - - Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb - character(len=20) :: name='csget' - logical :: rscale_, cscale_ - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - info = 0 - - nzin = 0 - if (present(imin)) then - imin_ = imin - else - imin_ = 1 - end if - if (present(imax)) then - imax_ = imax - else - imax_ = a%get_nrows() - end if - if (present(jmin)) then - jmin_ = jmin - else - jmin_ = 1 - end if - if (present(jmax)) then - jmax_ = jmax - else - jmax_ = a%get_ncols() - end if - if (present(rscale)) then - rscale_ = rscale - else - rscale_ = .true. - end if - if (present(cscale)) then - cscale_ = cscale - else - cscale_ = .true. - end if - - if (rscale_) then - mb = imax_ - imin_ +1 - else - mb = a%get_nrows() ! Should this be imax_ ?? - endif - if (cscale_) then - nb = jmax_ - jmin_ +1 - else - nb = a%get_ncols() ! Should this be jmax_ ?? - endif - call b%allocate(mb,nb) - - call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin_, jmax=jmax_, append=.false., & - & nzin=nzin, rscale=rscale_, cscale=cscale_) - - if (info /= 0) goto 9999 - - call b%set_nzeros(nzin+nzout) - call b%fix(info) - - if (info /= 0) goto 9999 - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_coo_csclip - - - !!$ !!$ subroutine d_coo_csget(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin) !!$ ! Output is always in COO format diff --git a/base/newserial/psbn_d_csr_mat_mod.f03 b/base/newserial/psbn_d_csr_mat_mod.f03 index 663b8ae9..c3a179e5 100644 --- a/base/newserial/psbn_d_csr_mat_mod.f03 +++ b/base/newserial/psbn_d_csr_mat_mod.f03 @@ -30,8 +30,6 @@ module psbn_d_csr_mat_mod procedure, pass(a) :: mv_to_fmt => d_mv_csr_to_fmt procedure, pass(a) :: mv_from_fmt => d_mv_csr_from_fmt procedure, pass(a) :: d_csgetrow => d_csr_csgetrow - procedure, pass(a) :: d_csgetblk => d_csr_csgetblk - procedure, pass(a) :: csclip => d_csr_csclip procedure, pass(a) :: get_size => d_csr_get_size procedure, pass(a) :: free => d_csr_free procedure, pass(a) :: trim => d_csr_trim @@ -44,8 +42,7 @@ module psbn_d_csr_mat_mod & d_mv_csr_to_coo, d_mv_csr_from_coo, & & d_cp_csr_to_fmt, d_cp_csr_from_fmt, & & d_mv_csr_to_fmt, d_mv_csr_from_fmt, & - & d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_csgetblk, & - & d_csr_csclip, d_csr_get_size + & d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_get_size interface diff --git a/prec/psb_ddiagsc_bld.f90 b/prec/psb_ddiagsc_bld.f90 index b7cdc49b..1b2ac00a 100644 --- a/prec/psb_ddiagsc_bld.f90 +++ b/prec/psb_ddiagsc_bld.f90 @@ -32,10 +32,11 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) use psb_base_mod + use psbn_d_mat_mod use psb_prec_mod, psb_protect_name => psb_ddiagsc_bld Implicit None - type(psb_dspmat_type), intent(in), target :: a + type(psbn_d_sparse_mat), intent(in), target :: a type(psb_desc_type), intent(in) :: desc_a type(psb_dprec_type),intent(inout) :: p character, intent(in) :: upd @@ -76,7 +77,7 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) ! ! Retrieve the diagonal entries of the matrix A ! - call psb_sp_getdiag(a,p%d,info) + call a%get_diag(p%d,info) if(info /= 0) then info=4010 ch_err='psb_sp_getdiag' @@ -105,18 +106,18 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) endif end do - if (a%pl(1) /= 0) then - ! - ! Apply the same row permutation as in the sparse matrix A - ! - call psb_gelp('n',a%pl,p%d,info) - if(info /= 0) then - info=4010 - ch_err='psb_dgelp' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif +!!$ if (a%pl(1) /= 0) then +!!$ ! +!!$ ! Apply the same row permutation as in the sparse matrix A +!!$ ! +!!$ call psb_gelp('n',a%pl,p%d,info) +!!$ if(info /= 0) then +!!$ info=4010 +!!$ ch_err='psb_dgelp' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ endif call psb_erractionrestore(err_act) return diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index 7962d1b3..68c5b711 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -102,7 +102,7 @@ subroutine psb_dprecbld(aa,desc_a,p,info,upd) case (psb_diag_) - call psb_diagsc_bld(a,desc_a,p,upd_,info) + call psb_diagsc_bld(aa,desc_a,p,upd_,info) if(info /= 0) then info=4010 ch_err='psb_diagsc_bld' diff --git a/prec/psb_dprecinit.f90 b/prec/psb_dprecinit.f90 index c79b52e3..513013aa 100644 --- a/prec/psb_dprecinit.f90 +++ b/prec/psb_dprecinit.f90 @@ -43,12 +43,7 @@ subroutine psb_dprecinit(p,ptype,info) call psb_realloc(psb_ifpsz,p%iprcparm,info) if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info) if (info /= 0) return - p%iprcparm(:) = 0 - p%iprcparm(:) = 0 - p%iprcparm(psb_p_type_) = psb_noprec_ - p%iprcparm(psb_f_type_) = psb_f_none_ - - return + select case(psb_toupper(ptype(1:len_trim(ptype)))) case ('NONE','NOPREC') p%iprcparm(:) = 0 @@ -60,12 +55,12 @@ subroutine psb_dprecinit(p,ptype,info) p%iprcparm(psb_p_type_) = psb_diag_ p%iprcparm(psb_f_type_) = psb_f_none_ - case ('BJAC') - p%iprcparm(:) = 0 - p%iprcparm(psb_p_type_) = psb_bjac_ - p%iprcparm(psb_f_type_) = psb_f_ilu_n_ - p%iprcparm(psb_ilu_fill_in_) = 0 - +!!$ case ('BJAC') +!!$ p%iprcparm(:) = 0 +!!$ p%iprcparm(psb_p_type_) = psb_bjac_ +!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_ +!!$ p%iprcparm(psb_ilu_fill_in_) = 0 +!!$ case default write(0,*) 'Unknown preconditioner type request "',ptype,'"' info = 2 diff --git a/prec/psb_prec_mod.f90 b/prec/psb_prec_mod.f90 index 67501ca4..d4e22b4d 100644 --- a/prec/psb_prec_mod.f90 +++ b/prec/psb_prec_mod.f90 @@ -406,8 +406,9 @@ module psb_prec_mod subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ use psb_prec_type, only : psb_dprec_type + use psbn_d_mat_mod integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a + type(psbn_d_sparse_mat), intent(in), target :: a type(psb_dprec_type), intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index e4e26360..0848bcf1 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -1,11 +1,11 @@ 7 Number of entries below this BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES -NONE Preconditioner NONE DIAG BJAC +DIAG Preconditioner NONE DIAG BJAC CSR Storage format for matrix A: CSR COO JAD -050 Domain size (acutal system is this**3) +060 Domain size (acutal system is this**3) 1 Stopping criterion -1000 MAXIT -040 ITRACE +0400 MAXIT +-40 ITRACE 20 IRST restart for RGMRES and BiCGSTABL diff --git a/test/serial/psbn_d_cxx_impl.f03 b/test/serial/psbn_d_cxx_impl.f03 index 1bc1b7c0..377202ed 100644 --- a/test/serial/psbn_d_cxx_impl.f03 +++ b/test/serial/psbn_d_cxx_impl.f03 @@ -1,16 +1,16 @@ -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! - ! - ! Computational routines - ! - ! - ! - ! - ! - ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!===================================== +! +! +! +! Computational routines +! +! +! +! +! +! +!===================================== subroutine d_cxx_csmv_impl(alpha,a,x,beta,y,info,trans) use psb_error_mod @@ -980,21 +980,199 @@ function d_cxx_csnmi_impl(a) result(res) end function d_cxx_csnmi_impl -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! - ! - ! Data management - ! - ! - ! - ! - ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!===================================== +! +! +! +! Data management +! +! +! +! +! +!===================================== + + +subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psbn_d_base_mat_mod + use psbn_d_cxx_mat_mod, psb_protect_name => d_cxx_csgetrow_impl + implicit none + + class(psbn_d_cxx_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_dpk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + logical :: append_, rscale_, cscale_ + integer :: nzin_, jmin_, jmax_, err_act, i + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax d_cxx_csput_impl @@ -1058,7 +1236,7 @@ contains use psb_const_mod use psb_realloc_mod use psb_string_mod - use psb_serial_mod + use psb_sort_mod implicit none class(psbn_d_cxx_sparse_mat), intent(inout) :: a @@ -1596,3 +1774,4 @@ subroutine d_cp_cxx_from_fmt_impl(a,b,info) if (info == 0) call a%mv_from_coo(tmp,info) end select end subroutine d_cp_cxx_from_fmt_impl + diff --git a/test/serial/psbn_d_cxx_mat_mod.f03 b/test/serial/psbn_d_cxx_mat_mod.f03 index 89d42ed1..b71bffbe 100644 --- a/test/serial/psbn_d_cxx_mat_mod.f03 +++ b/test/serial/psbn_d_cxx_mat_mod.f03 @@ -10,10 +10,13 @@ module psbn_d_cxx_mat_mod contains procedure, pass(a) :: get_nzeros => d_cxx_get_nzeros procedure, pass(a) :: get_fmt => d_cxx_get_fmt + procedure, pass(a) :: get_diag => d_cxx_get_diag procedure, pass(a) :: d_base_csmm => d_cxx_csmm procedure, pass(a) :: d_base_csmv => d_cxx_csmv procedure, pass(a) :: d_base_cssm => d_cxx_cssm procedure, pass(a) :: d_base_cssv => d_cxx_cssv + procedure, pass(a) :: d_scals => d_cxx_scals + procedure, pass(a) :: d_scal => d_cxx_scal procedure, pass(a) :: csnmi => d_cxx_csnmi procedure, pass(a) :: reallocate_nz => d_cxx_reallocate_nz procedure, pass(a) :: csput => d_cxx_csput @@ -26,16 +29,20 @@ module psbn_d_cxx_mat_mod procedure, pass(a) :: mv_from_coo => d_mv_cxx_from_coo procedure, pass(a) :: mv_to_fmt => d_mv_cxx_to_fmt procedure, pass(a) :: mv_from_fmt => d_mv_cxx_from_fmt + procedure, pass(a) :: d_csgetrow => d_cxx_csgetrow + procedure, pass(a) :: get_size => d_cxx_get_size procedure, pass(a) :: free => d_cxx_free + procedure, pass(a) :: trim => d_cxx_trim procedure, pass(a) :: print => d_cxx_print end type psbn_d_cxx_sparse_mat private :: d_cxx_get_nzeros, d_cxx_csmm, d_cxx_csmv, d_cxx_cssm, d_cxx_cssv, & & d_cxx_csput, d_cxx_reallocate_nz, d_cxx_allocate_mnnz, & - & d_cxx_free, d_cxx_print, d_cxx_get_fmt, d_cxx_csnmi, & + & d_cxx_free, d_cxx_print, d_cxx_get_fmt, d_cxx_csnmi, get_diag, & & d_cp_cxx_to_coo, d_cp_cxx_from_coo, & & d_mv_cxx_to_coo, d_mv_cxx_from_coo, & & d_cp_cxx_to_fmt, d_cp_cxx_from_fmt, & - & d_mv_cxx_to_fmt, d_mv_cxx_from_fmt + & d_mv_cxx_to_fmt, d_mv_cxx_from_fmt, & + & d_cxx_scals, d_cxx_scal, d_cxx_trim, d_cxx_csgetrow, d_cxx_get_size interface @@ -140,6 +147,26 @@ module psbn_d_cxx_mat_mod end subroutine d_cxx_csput_impl end interface + interface + subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psbn_d_cxx_sparse_mat + implicit none + + class(psbn_d_cxx_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_dpk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + end subroutine d_cxx_csgetrow_impl + end interface + interface d_cxx_cssm_impl subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans) use psb_const_mod @@ -195,7 +222,7 @@ module psbn_d_cxx_mat_mod contains -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== ! ! ! @@ -205,7 +232,7 @@ contains ! ! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== function d_cxx_get_fmt(a) result(res) implicit none @@ -221,8 +248,31 @@ contains res = a%irp(a%m+1)-1 end function d_cxx_get_nzeros + function d_cxx_get_size(a) result(res) + implicit none + class(psbn_d_cxx_sparse_mat), intent(in) :: a + integer :: res + + res = -1 + + if (allocated(a%ja)) then + if (res >= 0) then + res = min(res,size(a%ja)) + else + res = size(a%ja) + end if + end if + if (allocated(a%val)) then + if (res >= 0) then + res = min(res,size(a%val)) + else + res = size(a%val) + end if + end if + + end function d_cxx_get_size -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== ! ! ! @@ -232,7 +282,7 @@ contains ! ! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== subroutine d_cxx_reallocate_nz(nz,a) @@ -269,7 +319,7 @@ contains end subroutine d_cxx_reallocate_nz - subroutine d_cxx_csput(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + subroutine d_cxx_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) use psb_const_mod use psb_error_mod implicit none @@ -316,7 +366,7 @@ contains if (nz == 0) return - call d_cxx_csput_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + call d_cxx_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -332,6 +382,198 @@ contains return end subroutine d_cxx_csput + subroutine d_cxx_csgetrow(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psbn_d_cxx_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_dpk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + call d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_csgetrow + + + subroutine d_cxx_csgetblk(imin,imax,a,b,info,& + & jmin,jmax,iren,append,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psbn_d_cxx_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(inout) :: b + integer, intent(in) :: imin,imax + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax + logical, intent(in), optional :: rscale,cscale + Integer :: err_act, nzin, nzout + character(len=20) :: name='csget' + logical :: append_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(append)) then + append_ = append + else + append_ = .false. + endif + if (append_) then + nzin = a%get_nzeros() + else + nzin = 0 + endif + + call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin, jmax=jmax, iren=iren, append=append_, & + & nzin=nzin, rscale=rscale, cscale=cscale) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_csgetblk + + + subroutine d_cxx_csclip(a,b,info,& + & imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psbn_d_cxx_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer,intent(out) :: info + integer, intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb + character(len=20) :: name='csget' + logical :: rscale_, cscale_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + nzin = 0 + if (present(imin)) then + imin_ = imin + else + imin_ = 1 + end if + if (present(imax)) then + imax_ = imax + else + imax_ = a%get_nrows() + end if + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + end if + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + end if + if (present(rscale)) then + rscale_ = rscale + else + rscale_ = .true. + end if + if (present(cscale)) then + cscale_ = cscale + else + cscale_ = .true. + end if + + if (rscale_) then + mb = imax_ - imin_ +1 + else + mb = a%get_nrows() ! Should this be imax_ ?? + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = a%get_ncols() ! Should this be jmax_ ?? + endif + call b%allocate(mb,nb) + + call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=jmax_, append=.false., & + & nzin=nzin, rscale=rscale_, cscale=cscale_) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + + if (info /= 0) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_csclip + subroutine d_cxx_free(a) implicit none @@ -349,6 +591,38 @@ contains end subroutine d_cxx_free + subroutine d_cxx_trim(a) + use psb_realloc_mod + use psb_error_mod + implicit none + class(psbn_d_cxx_sparse_mat), intent(inout) :: a + Integer :: err_act, info, nz, m + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + m = a%get_nrows() + nz = a%get_nzeros() + if (info == 0) call psb_realloc(m+1,a%irp,info) + if (info == 0) call psb_realloc(nz,a%ja,info) + if (info == 0) call psb_realloc(nz,a%val,info) + + if (info /= 0) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_trim + subroutine d_cp_cxx_to_coo(a,b,info) use psb_error_mod @@ -671,7 +945,6 @@ contains subroutine d_cxx_print(iout,a,iv,eirs,eics,head,ivr,ivc) - use psb_spmat_type use psb_string_mod implicit none @@ -752,7 +1025,7 @@ contains end subroutine d_cxx_print -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== ! ! ! @@ -763,7 +1036,7 @@ contains ! ! ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !===================================== subroutine d_cxx_csmv(alpha,a,x,beta,y,info,trans) @@ -965,13 +1238,136 @@ contains character(len=20) :: name='csnmi' logical, parameter :: debug=.false. - + res = d_cxx_csnmi_impl(a) return end function d_cxx_csnmi + subroutine d_cxx_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_cxx_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act, mnm, i, j, k + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + + + do i=1, mnm + do k=a%irp(i),a%irp(i+1)-1 + j=a%ja(k) + if ((j==i) .and.(j <= mnm )) then + d(i) = a%val(k) + endif + enddo + end do + do i=mnm+1,size(d) + d(i) = dzero + end do + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_get_diag + + + subroutine d_cxx_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_cxx_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + m = a%get_nrows() + if (size(d) < m) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + + do i=1, m + do j = a%irp(i), a%irp(i+1) -1 + a%val(j) = a%val(j) * d(i) + end do + enddo + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_scal + + + subroutine d_cxx_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_cxx_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + + do i=1,a%get_nzeros() + a%val(i) = a%val(i) * d + enddo + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_cxx_scals + end module psbn_d_cxx_mat_mod diff --git a/util/psb_mat_dist_mod.f90 b/util/psb_mat_dist_mod.f90 index 6c6b26f7..fe5e75ef 100644 --- a/util/psb_mat_dist_mod.f90 +++ b/util/psb_mat_dist_mod.f90 @@ -547,13 +547,14 @@ contains ! on exit : unchanged. ! use psb_base_mod + use psbn_d_mat_mod implicit none ! parameters type(psb_dspmat_type) :: a_glob real(psb_dpk_) :: b_glob(:) integer :: ictxt - type(psb_dspmat_type) :: a + type(psbn_d_sparse_mat) :: a real(psb_dpk_), allocatable :: b(:) type(psb_desc_type) :: desc_a integer, intent(out) :: info