From fb28c925dc7a2dc2c49f2f2d94f72b0299662e6a Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 11 Sep 2009 14:20:06 +0000 Subject: [PATCH] psblas3: base/newserial/psbn_base_mat_mod.f03 base/newserial/psbn_d_base_mat_mod.f03 base/newserial/psbn_d_coo_impl.f03 base/newserial/psbn_d_csr_impl.f03 base/newserial/psbn_d_csr_mat_mod.f03 base/newserial/psbn_mat_mod.f03 test/serial/d_matgen.f03 Added SCAL. Defined CSGET/CSCLIP: they work on the same inner kernel implementation! --- base/newserial/psbn_base_mat_mod.f03 | 38 +- base/newserial/psbn_d_base_mat_mod.f03 | 620 ++++++++++++++++++++++++- base/newserial/psbn_d_coo_impl.f03 | 306 +++++++++++- base/newserial/psbn_d_csr_impl.f03 | 224 ++++++++- base/newserial/psbn_d_csr_mat_mod.f03 | 410 +++++++++++++++- base/newserial/psbn_mat_mod.f03 | 226 ++++++++- test/serial/d_matgen.f03 | 27 +- 7 files changed, 1774 insertions(+), 77 deletions(-) diff --git a/base/newserial/psbn_base_mat_mod.f03 b/base/newserial/psbn_base_mat_mod.f03 index 04c5a964..a8fa96d0 100644 --- a/base/newserial/psbn_base_mat_mod.f03 +++ b/base/newserial/psbn_base_mat_mod.f03 @@ -84,6 +84,7 @@ module psbn_base_mat_mod procedure, pass(a) :: allocate_mnnz procedure, pass(a) :: reallocate_nz procedure, pass(a) :: free + procedure, pass(a) :: trim generic, public :: allocate => allocate_mnnz generic, public :: reallocate => reallocate_nz procedure, pass(a) :: print => sparse_print @@ -96,7 +97,7 @@ module psbn_base_mat_mod & get_nzeros, get_size, get_state, get_dupl, is_null, is_bld, & & is_upd, is_asb, is_sorted, is_upper, is_lower, is_triangle, & & is_unit, get_neigh, allocate_mn, allocate_mnnz, reallocate_nz, & - & free, sparse_print,get_fmt + & free, sparse_print, get_fmt, trim contains @@ -329,7 +330,7 @@ contains character(len=20) :: name='base_get_nzeros' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + call psb_get_erraction(err_act) res = -1 ! This is the base version. If we get here ! it means the derived class is incomplete, @@ -353,7 +354,7 @@ contains character(len=20) :: name='get_size' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + call psb_get_erraction(err_act) res = -1 ! This is the base version. If we get here ! it means the derived class is incomplete, @@ -383,7 +384,7 @@ contains character(len=20) :: name='sparse_print' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + call psb_get_erraction(err_act) info = 700 ! This is the base version. If we get here ! it means the derived class is incomplete, @@ -412,7 +413,7 @@ contains character(len=20) :: name='get_neigh' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + call psb_get_erraction(err_act) info = 700 ! This is the base version. If we get here ! it means the derived class is incomplete, @@ -436,7 +437,7 @@ contains character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + 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. @@ -458,7 +459,7 @@ contains character(len=20) :: name='reallocate_nz' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + 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. @@ -479,7 +480,7 @@ contains character(len=20) :: name='free' logical, parameter :: debug=.false. - call psb_erractionsave(err_act) + 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. @@ -492,5 +493,26 @@ contains end subroutine free + subroutine trim(a) + use psb_error_mod + implicit none + class(psbn_base_sparse_mat), intent(inout) :: a + Integer :: err_act + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + 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. + call psb_errpush(700,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine trim + end module psbn_base_mat_mod diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 8f51d99f..db43d352 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -10,8 +10,16 @@ module psbn_d_base_mat_mod procedure, pass(a) :: d_base_cssv procedure, pass(a) :: d_base_cssm generic, public :: cssm => d_base_cssm, d_base_cssv + procedure, pass(a) :: d_scals + procedure, pass(a) :: d_scal + generic, public :: scal => d_scals, d_scal + procedure, pass(a) :: get_diag procedure, pass(a) :: csnmi procedure, pass(a) :: csput + procedure, pass(a) :: d_csgetrow + procedure, pass(a) :: d_csgetblk + generic, public :: csget => d_csgetrow, d_csgetblk + procedure, pass(a) :: csclip procedure, pass(a) :: cp_to_coo procedure, pass(a) :: cp_from_coo procedure, pass(a) :: cp_to_fmt @@ -21,9 +29,12 @@ module psbn_d_base_mat_mod procedure, pass(a) :: mv_to_fmt procedure, pass(a) :: mv_from_fmt end type psbn_d_base_sparse_mat + private :: d_base_csmv, d_base_csmm, d_base_cssv, d_base_cssm,& - & csnmi, csput, cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, & - & mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt + & d_scals, d_scal, csnmi, csput, d_csgetrow, d_csgetblk, & + & cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, & + & mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt, & + & get_diag, csclip type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat @@ -41,8 +52,11 @@ module psbn_d_base_mat_mod procedure, pass(a) :: d_base_csmv => d_coo_csmv procedure, pass(a) :: d_base_cssm => d_coo_cssm procedure, pass(a) :: d_base_cssv => d_coo_cssv + procedure, pass(a) :: d_scals => d_coo_scals + procedure, pass(a) :: d_scal => d_coo_scal procedure, pass(a) :: csnmi => d_coo_csnmi procedure, pass(a) :: csput => d_coo_csput + procedure, pass(a) :: get_diag => d_coo_get_diag procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz procedure, pass(a) :: allocate_mnnz => d_coo_allocate_mnnz procedure, pass(a) :: cp_to_coo => d_cp_coo_to_coo @@ -55,16 +69,23 @@ module psbn_d_base_mat_mod procedure, pass(a) :: mv_from_fmt => d_mv_coo_from_fmt procedure, pass(a) :: fix => d_fix_coo 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 end type psbn_d_coo_sparse_mat - private :: d_coo_get_nzeros, d_coo_set_nzeros, & + + private :: d_coo_get_nzeros, d_coo_set_nzeros, d_coo_get_diag, & & d_coo_csmm, d_coo_csmv, d_coo_cssm, d_coo_cssv, d_coo_csnmi, & & d_coo_csput, d_coo_reallocate_nz, d_coo_allocate_mnnz, & & 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_cp_coo_to_fmt, d_cp_coo_from_fmt, & + & d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_csgetblk, & + & d_coo_csclip interface @@ -175,12 +196,33 @@ module psbn_d_base_mat_mod import psbn_d_coo_sparse_mat class(psbn_d_coo_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: val(:) - integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(in) :: nz,ia(:), ja(:),& + & imin,imax,jmin,jmax integer, intent(out) :: info integer, intent(in), optional :: gtl(:) end subroutine d_coo_csput_impl end interface + interface + subroutine d_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psbn_d_coo_sparse_mat + implicit none + + class(psbn_d_coo_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_coo_csgetrow_impl + end interface + interface d_coo_cssm_impl subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans) use psb_const_mod @@ -490,6 +532,108 @@ contains end subroutine csput + subroutine d_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_base_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_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 (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine d_csgetrow + + + subroutine d_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_base_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 + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + 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 (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine d_csgetblk + + + subroutine 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_base_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 + character(len=20) :: name='csclip' + logical, parameter :: debug=.false. + + 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 (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine csclip + + + !==================================== ! ! @@ -610,9 +754,61 @@ contains end if return - end subroutine d_base_cssv + + subroutine d_scals(d,a,info) + use psb_error_mod + implicit none + class(psbn_d_base_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='d_scals' + logical, parameter :: debug=.false. + + 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 (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine d_scals + + + subroutine d_scal(d,a,info) + use psb_error_mod + implicit none + class(psbn_d_base_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='d_scal' + logical, parameter :: debug=.false. + + 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 (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine d_scal + + function csnmi(a) result(res) use psb_error_mod use psb_const_mod @@ -640,6 +836,33 @@ contains end function csnmi + subroutine get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_base_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + 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 (err_act /= psb_act_ret_) then + call psb_error() + end if + + return + + end subroutine get_diag + @@ -1125,6 +1348,242 @@ contains end subroutine d_coo_csput + subroutine d_coo_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_coo_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_coo_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_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 +!!$ use psb_error_mod +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ class(psbn_d_coo_sparse_mat), intent(inout) :: a +!!$ integer, intent(in) :: irw +!!$ 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 :: lrw, nzin +!!$ Integer :: err_act +!!$ character(len=20) :: name='csget' +!!$ logical, parameter :: debug=.false. +!!$ +!!$ call psb_erractionsave(err_act) +!!$ info = 0 +!!$ +!!$ call d_coo_csget_impl(irw,a,nz,ia,ja,val,info,iren,lrw,append,nzin) +!!$ 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_csget +!!$ + + subroutine d_coo_free(a) implicit none @@ -1141,6 +1600,38 @@ contains end subroutine d_coo_free + + subroutine d_coo_trim(a) + use psb_realloc_mod + use psb_error_mod + implicit none + class(psbn_d_coo_sparse_mat), intent(inout) :: a + Integer :: err_act, info, nz + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + nz = a%get_nzeros() + if (info == 0) call psb_realloc(nz,a%ia,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_coo_trim + subroutine d_coo_allocate_mnnz(m,n,a,nz) use psb_error_mod use psb_realloc_mod @@ -1275,6 +1766,7 @@ contains + !==================================== ! ! @@ -1494,6 +1986,122 @@ contains end function d_coo_csnmi + subroutine d_coo_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j + 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 + d(:) = dzero + + do i=1,a%get_nzeros() + j=a%ia(i) + if ((j==a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then + d(j) = a%val(i) + endif + 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_coo_get_diag + + subroutine d_coo_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_coo_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,a%get_nzeros() + j = a%ia(i) + a%val(i) = a%val(i) * d(j) + 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_coo_scal + + subroutine d_coo_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_coo_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_coo_scals + end module psbn_d_base_mat_mod diff --git a/base/newserial/psbn_d_coo_impl.f03 b/base/newserial/psbn_d_coo_impl.f03 index de79522e..58ac106a 100644 --- a/base/newserial/psbn_d_coo_impl.f03 +++ b/base/newserial/psbn_d_coo_impl.f03 @@ -881,26 +881,309 @@ function d_coo_csnmi_impl(a) result(res) acc = acc + abs(a%val(k)) end do res = max(res,acc) + i = j end do end function d_coo_csnmi_impl -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! - ! - ! Data management - ! - ! - ! - ! - ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!==================================== +! +! +! +! Data management +! +! +! +! +! +!==================================== +subroutine d_coo_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, psb_protect_name => d_coo_csgetrow_impl + implicit none + + class(psbn_d_coo_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= psb_debug_serial_)& + & write(debug_unit,*) trim(name), ': srtdcoo ' + do + ip = psb_ibsrch(irw,nza,a%ia) + if (ip /= -1) exit + irw = irw + 1 + if (irw > imax) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error? ',& + & irw,lrw,imin + exit + end if + end do + + if (ip /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (ip < 2) exit + if (a%ia(ip-1) == irw) then + ip = ip -1 + else + exit + end if + end do + + end if + + do + jp = psb_ibsrch(lrw,nza,a%ia) + if (jp /= -1) exit + lrw = lrw - 1 + if (irw > lrw) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error?' + exit + end if + end do + + if (jp /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (jp == nza) exit + if (a%ia(jp+1) == lrw) then + jp = jp + 1 + else + exit + end if + end do + end if + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': ip jp',ip,jp,nza + if ((ip /= -1) .and.(jp /= -1)) then + ! Now do the copy. + nzt = jp - ip +1 + nz = 0 + + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + if (present(iren)) then + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + val(nzin_) = a%val(i) + ia(nzin_) = iren(a%ia(i)) + ja(nzin_) = iren(a%ja(i)) + end if + enddo + else + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + val(nzin_) = a%val(i) + ia(nzin_) = a%ia(i) + ja(nzin_) = a%ja(i) + end if + enddo + end if + else + nz = 0 + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': unsorted ' + + nzt = (nza*(lrw-irw+1))/max(a%get_nrows(),1) + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + if (present(iren)) then + k = 0 + do i=1, a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + end if + val(nzin_+k) = a%val(i) + ia(nzin_+k) = iren(a%ia(i)) + ja(nzin_+k) = iren(a%ja(i)) + endif + enddo + else + k = 0 + do i=1,a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + end if + val(nzin_+k) = a%val(i) + ia(nzin_+k) = (a%ia(i)) + ja(nzin_+k) = (a%ja(i)) + endif + enddo + nzin_=nzin_+k + end if + nz = k + end if + + end subroutine coo_getrow + +end subroutine d_coo_csgetrow_impl + + subroutine d_coo_csput_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) use psb_error_mod use psb_realloc_mod @@ -965,6 +1248,7 @@ subroutine d_coo_csput_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,& & imin,imax,jmin,jmax,info,gtl) call a%set_nzeros(nza) + call a%set_sorted(.false.) else if (a%is_upd()) then diff --git a/base/newserial/psbn_d_csr_impl.f03 b/base/newserial/psbn_d_csr_impl.f03 index fe0181a6..e9545baa 100644 --- a/base/newserial/psbn_d_csr_impl.f03 +++ b/base/newserial/psbn_d_csr_impl.f03 @@ -1,16 +1,16 @@ -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! - ! - ! Computational routines - ! - ! - ! - ! - ! - ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!===================================== +! +! +! +! Computational routines +! +! +! +! +! +! +!===================================== subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans) use psb_error_mod @@ -980,17 +980,195 @@ function d_csr_csnmi_impl(a) result(res) end function d_csr_csnmi_impl -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! - ! - ! Data management - ! - ! - ! - ! - ! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!===================================== +! +! +! +! Data management +! +! +! +! +! +!===================================== + + +subroutine d_csr_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_csr_mat_mod, psb_protect_name => d_csr_csgetrow_impl + implicit none + + class(psbn_d_csr_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_csr_get_nzeros procedure, pass(a) :: get_fmt => d_csr_get_fmt + procedure, pass(a) :: get_diag => d_csr_get_diag procedure, pass(a) :: d_base_csmm => d_csr_csmm procedure, pass(a) :: d_base_csmv => d_csr_csmv procedure, pass(a) :: d_base_cssm => d_csr_cssm procedure, pass(a) :: d_base_cssv => d_csr_cssv + procedure, pass(a) :: d_scals => d_csr_scals + procedure, pass(a) :: d_scal => d_csr_scal procedure, pass(a) :: csnmi => d_csr_csnmi procedure, pass(a) :: reallocate_nz => d_csr_reallocate_nz procedure, pass(a) :: csput => d_csr_csput @@ -26,16 +29,23 @@ module psbn_d_csr_mat_mod procedure, pass(a) :: mv_from_coo => d_mv_csr_from_coo 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 procedure, pass(a) :: print => d_csr_print end type psbn_d_csr_sparse_mat private :: d_csr_get_nzeros, d_csr_csmm, d_csr_csmv, d_csr_cssm, d_csr_cssv, & & d_csr_csput, d_csr_reallocate_nz, d_csr_allocate_mnnz, & - & d_csr_free, d_csr_print, d_csr_get_fmt, d_csr_csnmi, & + & d_csr_free, d_csr_print, d_csr_get_fmt, d_csr_csnmi, get_diag, & & d_cp_csr_to_coo, d_cp_csr_from_coo, & & 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_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 interface @@ -140,6 +150,26 @@ module psbn_d_csr_mat_mod end subroutine d_csr_csput_impl end interface + interface + subroutine d_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psbn_d_csr_sparse_mat + implicit none + + class(psbn_d_csr_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_csr_csgetrow_impl + end interface + interface d_csr_cssm_impl subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) use psb_const_mod @@ -221,8 +251,31 @@ contains res = a%irp(a%m+1)-1 end function d_csr_get_nzeros + function d_csr_get_size(a) result(res) + implicit none + class(psbn_d_csr_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_csr_get_size -!===================================== + !===================================== ! ! ! @@ -232,7 +285,7 @@ contains ! ! ! -!===================================== + !===================================== subroutine d_csr_reallocate_nz(nz,a) @@ -332,6 +385,198 @@ contains return end subroutine d_csr_csput + subroutine d_csr_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_csr_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_csr_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_csr_csgetrow + + + subroutine d_csr_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_csr_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_csr_csgetblk + + + subroutine d_csr_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_csr_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_csr_csclip + subroutine d_csr_free(a) implicit none @@ -349,6 +594,38 @@ contains end subroutine d_csr_free + subroutine d_csr_trim(a) + use psb_realloc_mod + use psb_error_mod + implicit none + class(psbn_d_csr_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_csr_trim + subroutine d_cp_csr_to_coo(a,b,info) use psb_error_mod @@ -965,13 +1242,136 @@ contains character(len=20) :: name='csnmi' logical, parameter :: debug=.false. - + res = d_csr_csnmi_impl(a) return end function d_csr_csnmi + subroutine d_csr_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_csr_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_csr_get_diag + + + subroutine d_csr_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_csr_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_csr_scal + + + subroutine d_csr_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_csr_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_csr_scals + end module psbn_d_csr_mat_mod diff --git a/base/newserial/psbn_mat_mod.f03 b/base/newserial/psbn_mat_mod.f03 index 9b9a47e3..81548c93 100644 --- a/base/newserial/psbn_mat_mod.f03 +++ b/base/newserial/psbn_mat_mod.f03 @@ -44,8 +44,12 @@ module psbn_d_mat_mod ! Memory/data management procedure, pass(a) :: csall procedure, pass(a) :: free + procedure, pass(a) :: trim procedure, pass(a) :: csput - procedure, pass(a) :: csget + procedure, pass(a) :: d_csgetrow + procedure, pass(a) :: d_csgetblk + generic, public :: csget => d_csgetrow, d_csgetblk + procedure, pass(a) :: csclip procedure, pass(a) :: reall => reallocate_nz procedure, pass(a) :: get_neigh procedure, pass(a) :: d_cscnv @@ -54,6 +58,7 @@ module psbn_d_mat_mod procedure, pass(a) :: print => sparse_print ! Computational routines + procedure, pass(a) :: get_diag procedure, pass(a) :: csnmi procedure, pass(a) :: d_csmv procedure, pass(a) :: d_csmm @@ -68,11 +73,12 @@ module psbn_d_mat_mod private :: get_nrows, get_ncols, get_nzeros, get_size, & & get_state, get_dupl, is_null, is_bld, is_upd, & & is_asb, is_sorted, is_upper, is_lower, is_triangle, & - & is_unit, get_neigh, csall, csput, csget, d_cscnv, d_cscnv_ip, & - & reallocate_nz, free, d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, & + & is_unit, get_neigh, csall, csput, d_csgetrow,& + & d_csgetblk, csclip, d_cscnv, d_cscnv_ip, & + & reallocate_nz, free, trim, d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, & & set_nrows, set_ncols, set_dupl, set_state, set_null, set_bld, & & set_upd, set_asb, set_sorted, set_upper, set_lower, set_triangle, & - & set_unit, csnmi + & set_unit, csnmi, get_diag contains @@ -904,6 +910,7 @@ contains character(len=20) :: name='reallocate_nz' logical, parameter :: debug=.false. + call psb_get_erraction(err_act) if (.not.allocated(a%a)) then info = 1121 call psb_errpush(info,name) @@ -933,6 +940,7 @@ contains character(len=20) :: name='free' logical, parameter :: debug=.false. + call psb_get_erraction(err_act) if (.not.allocated(a%a)) then info = 1121 call psb_errpush(info,name) @@ -944,7 +952,6 @@ contains return 9999 continue - call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then call psb_error() @@ -954,6 +961,36 @@ contains end subroutine free + subroutine trim(a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%trim() + + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine trim + + subroutine csput(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) use psbn_d_base_mat_mod use psb_error_mod @@ -993,20 +1030,79 @@ contains end subroutine csput - subroutine csget(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + subroutine d_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 use psbn_d_base_mat_mod + implicit none + + class(psbn_d_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. + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call a%a%csget(imin,imax,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 + + end subroutine d_csgetrow + + + + subroutine d_csgetblk(imin,imax,a,b,info,& + & jmin,jmax,iren,append,rscale,cscale) + ! Output is always in COO format use psb_error_mod - implicit none - class(psbn_d_sparse_mat), intent(inout) :: a - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: nz, ia(:), ja(:) - integer, intent(in) :: imin,imax,jmin,jmax - integer, intent(out) :: info - integer, intent(in), optional :: gtl(:) + use psb_const_mod + use psbn_d_base_mat_mod + implicit none + + class(psbn_d_sparse_mat), intent(in) :: a + class(psbn_d_sparse_mat), intent(out) :: 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 - character(len=20) :: name='csput' + character(len=20) :: name='csget' logical, parameter :: debug=.false. + type(psbn_d_coo_sparse_mat), allocatable :: acoo + info = 0 call psb_erractionsave(err_act) @@ -1016,14 +1112,60 @@ contains goto 9999 endif - info = 700 - call psb_errpush(info,name,a_err='CSGET') - goto 9999 + allocate(acoo,stat=info) + + if (info == 0) call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) + if (info == 0) call move_alloc(acoo,b%a) + 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 + + end subroutine d_csgetblk -!!$ -!!$ call a%a%csget(nz,val,ia,ja,imin,imax,jmin,jmax,info,gtl) -!!$ if (info /= 0) goto 9999 + + subroutine csclip(a,b,info,& + & imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psbn_d_base_mat_mod + implicit none + + class(psbn_d_sparse_mat), intent(in) :: a + class(psbn_d_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 + character(len=20) :: name='csclip' + logical, parameter :: debug=.false. + type(psbn_d_coo_sparse_mat), allocatable :: acoo + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + if (info == 0) call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) + if (info == 0) call move_alloc(acoo,b%a) + if (info /= 0) goto 9999 call psb_erractionrestore(err_act) return @@ -1036,7 +1178,8 @@ contains return end if - end subroutine csget + end subroutine csclip + subroutine d_cscnv(a,b,info,type,mold,upd,dupl) @@ -1399,6 +1542,7 @@ contains character(len=20) :: name='csnmi' logical, parameter :: debug=.false. + call psb_get_erraction(err_act) if (.not.allocated(a%a)) then info = 1121 call psb_errpush(info,name) @@ -1420,5 +1564,45 @@ contains end function csnmi + + + subroutine get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%get_diag(d,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 get_diag + + + end module psbn_d_mat_mod diff --git a/test/serial/d_matgen.f03 b/test/serial/d_matgen.f03 index e8f9b6de..2733fe41 100644 --- a/test/serial/d_matgen.f03 +++ b/test/serial/d_matgen.f03 @@ -1,5 +1,5 @@ ! -program d_coo_matgen +program d_matgen use psb_base_mod use psb_prec_mod use psb_krylov_mod @@ -154,7 +154,7 @@ contains integer :: np, iam, nr, nt,nz,isz integer :: element integer, allocatable :: irow(:),icol(:),myidx(:) - real(psb_dpk_), allocatable :: val(:) + real(psb_dpk_), allocatable :: val(:), diag(:) type(psbn_d_sparse_mat) :: a_n type(psbn_d_coo_sparse_mat) :: acoo type(psbn_d_csr_sparse_mat) :: acsr @@ -382,6 +382,26 @@ contains call a_n%print(20) anorm = a_n%csnmi() write(0,*) 'Nrm infinity ',anorm + call a_n%csget(2,3,element,irow,icol,val,info) + write(0,*) 'From csget ',element,info + if (info == 0) then + do i=1,element + write(0,*) irow(i),icol(i),val(i) + end do + end if + + isz = a_n%get_size() + write(0,*) 'Size 1: ',isz + call a_n%trim() + isz = a_n%get_size() + write(0,*) 'Size 2: ',isz + + + + allocate(diag(nlr),stat=info) + if (info == 0) then + call a_n%get_diag(diag,info) + end if !!$ t1 = psb_wtime() call a_n%cscnv(info,mold=acxx) @@ -396,6 +416,7 @@ contains call a_n%print(21) anorm = a_n%csnmi() write(0,*) 'Nrm infinity ',anorm + !!$ if(iam == psb_root_) then @@ -419,7 +440,7 @@ contains end if return end subroutine create_matrix -end program d_coo_matgen +end program d_matgen ! ! functions parametrizing the differential equation !