diff --git a/Changelog b/Changelog index 041653a3..b7f3cd21 100644 --- a/Changelog +++ b/Changelog @@ -1,6 +1,8 @@ Changelog. A lot less detailed than usual, at least for past history. - +2009/09/15: First working OO implementation for serial routines on sparse + matrix data structures. Only D for the time being. + 2009/08/25: New configure flag --enable-serial for serial-only compilation. diff --git a/base/modules/psb_base_mat_mod.f03 b/base/modules/psb_base_mat_mod.f03 index 7af04cad..32047159 100644 --- a/base/modules/psb_base_mat_mod.f03 +++ b/base/modules/psb_base_mat_mod.f03 @@ -2,30 +2,8 @@ module psb_base_mat_mod use psb_const_mod -!!$ integer, parameter :: psb_invalid_ = -1 -!!$ integer, parameter :: psb_spmat_null_=0, psb_spmat_bld_=1 -!!$ integer, parameter :: psb_spmat_asb_=2, psb_spmat_upd_=4 -!!$ -!!$ integer, parameter :: psb_ireg_flgs_=10, psb_ip2_=0 -!!$ integer, parameter :: psb_iflag_=2, psb_ichk_=3 -!!$ integer, parameter :: psb_nnzt_=4, psb_zero_=5,psb_ipc_=6 -!!$ ! Duplicate coefficients handling -!!$ ! These are usually set while calling spcnv as one of its -!!$ ! optional arugments. -!!$ integer, parameter :: psb_dupl_ovwrt_ = 0 -!!$ integer, parameter :: psb_dupl_add_ = 1 -!!$ integer, parameter :: psb_dupl_err_ = 2 -!!$ integer, parameter :: psb_dupl_def_ = psb_dupl_ovwrt_ -!!$ ! Matrix update mode -!!$ integer, parameter :: psb_upd_srch_ = 98764 -!!$ integer, parameter :: psb_upd_perm_ = 98765 -!!$ integer, parameter :: psb_upd_dflt_ = psb_upd_srch_ -!!$ integer, parameter :: psb_maxjdrows_=8, psb_minjdrows_=4 -!!$ integer, parameter :: psb_dbleint_=2 - - type :: psb_base_sparse_mat - integer :: m, n + integer, private :: m, n integer, private :: state, duplicate logical, private :: triangle, unitd, upper, sorted contains @@ -39,6 +17,7 @@ module psb_base_mat_mod procedure, pass(a) :: get_nrows procedure, pass(a) :: get_ncols procedure, pass(a) :: get_nzeros + procedure, pass(a) :: get_nz_row procedure, pass(a) :: get_size procedure, pass(a) :: get_state procedure, pass(a) :: get_dupl @@ -84,8 +63,11 @@ module psb_base_mat_mod procedure, pass(a) :: reallocate_nz procedure, pass(a) :: free procedure, pass(a) :: trim + procedure, pass(a) :: reinit generic, public :: allocate => allocate_mnnz generic, public :: reallocate => reallocate_nz + procedure, pass(a) :: csgetptn + generic, public :: csget => csgetptn procedure, pass(a) :: print => sparse_print procedure, pass(a) :: sizeof @@ -97,12 +79,11 @@ module psb_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, trim, sizeof + & free, sparse_print, get_fmt, trim, sizeof, reinit, csgetptn, & + & get_nz_row - contains - function sizeof(a) result(res) implicit none class(psb_base_sparse_mat), intent(in) :: a @@ -331,6 +312,31 @@ contains end function is_sorted + function get_nz_row(idx,a) result(res) + use psb_error_mod + implicit none + integer, intent(in) :: idx + class(psb_base_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act + character(len=20) :: name='base_get_nz_row' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + res = -1 + ! 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 function get_nz_row + function get_nzeros(a) result(res) use psb_error_mod implicit none @@ -379,6 +385,76 @@ contains end function get_size + subroutine reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + + Integer :: err_act, info + character(len=20) :: name='reinit' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + info = 700 + ! 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 reinit + +!!$ +!!$ ! +!!$ ! Since at this level we have only simple components, +!!$ ! mv_from is identical to cp_from. +!!$ ! +!!$ subroutine mv_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_base_sparse_mat), intent(out) :: a +!!$ type(psb_base_sparse_mat), intent(inout) :: b +!!$ +!!$ a%m = b%m +!!$ a%n = b%n +!!$ a%state = b%state +!!$ a%duplicate = b%duplicate +!!$ a%triangle = b%triangle +!!$ a%unitd = b%unitd +!!$ a%upper = b%upper +!!$ a%sorted = b%sorted +!!$ +!!$ return +!!$ +!!$ end subroutine mv_from +!!$ +!!$ subroutine cp_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_base_sparse_mat), intent(out) :: a +!!$ type(psb_base_sparse_mat), intent(in) :: b +!!$ +!!$ a%m = b%m +!!$ a%n = b%n +!!$ a%state = b%state +!!$ a%duplicate = b%duplicate +!!$ a%triangle = b%triangle +!!$ a%unitd = b%unitd +!!$ a%upper = b%upper +!!$ a%sorted = b%sorted +!!$ +!!$ return +!!$ +!!$ end subroutine cp_from +!!$ subroutine sparse_print(iout,a,iv,eirs,eics,head,ivr,ivc) use psb_error_mod implicit none @@ -408,6 +484,39 @@ contains end subroutine sparse_print + subroutine csgetptn(imin,imax,a,nz,ia,ja,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(psb_base_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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 csgetptn subroutine get_neigh(a,idx,neigh,n,info,lev) use psb_error_mod @@ -524,5 +633,6 @@ contains end subroutine trim + end module psb_base_mat_mod diff --git a/base/modules/psb_base_mod.f90 b/base/modules/psb_base_mod.f90 index c6e58a8a..03faccc5 100644 --- a/base/modules/psb_base_mod.f90 +++ b/base/modules/psb_base_mod.f90 @@ -41,4 +41,5 @@ module psb_base_mod use psb_psblas_mod use psb_gps_mod use psb_tools_mod + use psb_mat_mod end module psb_base_mod diff --git a/base/modules/psb_d_base_mat_mod.f03 b/base/modules/psb_d_base_mat_mod.f03 index f4a7baef..71dc4f3d 100644 --- a/base/modules/psb_d_base_mat_mod.f03 +++ b/base/modules/psb_d_base_mat_mod.f03 @@ -38,9 +38,8 @@ module psb_d_base_mat_mod & 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, d_cssv, d_cssm - - + & get_diag, csclip, d_cssv, d_cssm + type, extends(psb_d_base_sparse_mat) :: psb_d_coo_sparse_mat integer :: nnz @@ -75,9 +74,12 @@ module psb_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) :: csgetptn => d_coo_csgetptn procedure, pass(a) :: print => d_coo_print procedure, pass(a) :: get_fmt => d_coo_get_fmt + procedure, pass(a) :: get_nz_row => d_coo_get_nz_row procedure, pass(a) :: sizeof => d_coo_sizeof + procedure, pass(a) :: reinit => d_coo_reinit end type psb_d_coo_sparse_mat @@ -87,8 +89,9 @@ module psb_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_sizeof - + & d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_sizeof, & + & d_coo_csgetptn, d_coo_get_nz_row, d_coo_reinit + interface subroutine d_fix_coo_inner(nzin,dupl,ia,ja,val,nzout,info,idir) @@ -125,7 +128,7 @@ module psb_d_base_mat_mod subroutine d_cp_coo_from_coo_impl(a,b,info) use psb_const_mod import psb_d_coo_sparse_mat - class(psb_d_coo_sparse_mat), intent(inout) :: a + class(psb_d_coo_sparse_mat), intent(out) :: a class(psb_d_coo_sparse_mat), intent(in) :: b integer, intent(out) :: info end subroutine d_cp_coo_from_coo_impl @@ -205,6 +208,24 @@ module psb_d_base_mat_mod end subroutine d_coo_csput_impl end interface + interface + subroutine d_coo_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_d_coo_sparse_mat + implicit none + class(psb_d_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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_csgetptn_impl + end interface + interface subroutine d_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) @@ -293,7 +314,6 @@ contains ! !==================================== - subroutine cp_to_coo(a,b,info) use psb_error_mod use psb_realloc_mod @@ -1278,6 +1298,56 @@ contains end function d_coo_get_nzeros + function d_coo_get_nz_row(idx,a) result(res) + use psb_const_mod + use psb_sort_mod + implicit none + + class(psb_d_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer :: res + integer :: nzin_, nza,ip,jp,i,k + + res = 0 + nza = a%get_nzeros() + if (a%is_sorted()) then + ! In this case we can do a binary search. + ip = psb_ibsrch(idx,nza,a%ia) + if (ip /= -1) return + jp = ip + do + if (ip < 2) exit + if (a%ia(ip-1) == idx) then + ip = ip -1 + else + exit + end if + end do + do + if (jp == nza) exit + if (a%ia(jp+1) == idx) then + jp = jp + 1 + else + exit + end if + end do + + res = jp - ip +1 + + else + + res = 0 + + do i=1, nza + if (a%ia(i) == idx) then + res = res + 1 + end if + end do + + end if + + end function d_coo_get_nz_row + !==================================== ! ! @@ -1383,7 +1453,7 @@ contains use psb_error_mod use psb_realloc_mod implicit none - class(psb_d_coo_sparse_mat), intent(inout) :: a + class(psb_d_coo_sparse_mat), intent(out) :: a class(psb_d_coo_sparse_mat), intent(in) :: b integer, intent(out) :: info @@ -1411,6 +1481,7 @@ contains end subroutine d_cp_coo_from_coo + subroutine d_cp_coo_to_fmt(a,b,info) use psb_error_mod use psb_realloc_mod @@ -1541,6 +1612,73 @@ contains end subroutine d_mv_coo_from_coo + +!!$ +!!$ subroutine d_coo_cp_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_d_coo_sparse_mat), intent(out) :: a +!!$ type(psb_d_coo_sparse_mat), intent(inout) :: b +!!$ +!!$ +!!$ Integer :: err_act, info +!!$ character(len=20) :: name='mv_from' +!!$ logical, parameter :: debug=.false. +!!$ +!!$ call psb_erractionsave(err_act) +!!$ info = 0 +!!$ call d_cp_coo_from_coo_impl(a,b,info) +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 continue +!!$ call psb_erractionrestore(err_act) +!!$ +!!$ call psb_errpush(info,name) +!!$ +!!$ if (err_act /= psb_act_ret_) then +!!$ call psb_error() +!!$ end if +!!$ return +!!$ +!!$ end subroutine d_coo_cp_from +!!$ +!!$ subroutine d_coo_mv_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_d_coo_sparse_mat), intent(out) :: a +!!$ type(psb_d_coo_sparse_mat), intent(inout) :: b +!!$ +!!$ +!!$ Integer :: err_act, info +!!$ character(len=20) :: name='mv_from' +!!$ logical, parameter :: debug=.false. +!!$ +!!$ call psb_erractionsave(err_act) +!!$ info = 0 +!!$ call d_mv_coo_from_coo_impl(a,b,info) +!!$ if (info /= 0) goto 9999 +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 continue +!!$ call psb_erractionrestore(err_act) +!!$ +!!$ call psb_errpush(info,name) +!!$ +!!$ if (err_act /= psb_act_ret_) then +!!$ call psb_error() +!!$ end if +!!$ return +!!$ +!!$ end subroutine d_coo_mv_from +!!$ + subroutine d_mv_coo_to_fmt(a,b,info) use psb_error_mod use psb_realloc_mod @@ -1750,46 +1888,47 @@ contains end subroutine d_coo_csgetrow -!!$ -!!$ 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(psb_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_csgetptn(imin,imax,a,nz,ia,ja,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(psb_d_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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_csgetptn_impl(imin,imax,a,nz,ia,ja,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_csgetptn subroutine d_coo_free(a) @@ -1808,6 +1947,54 @@ contains end subroutine d_coo_free + subroutine d_coo_reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_d_coo_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + + Integer :: err_act, info + character(len=20) :: name='reinit' + logical :: clear_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + + if (present(clear)) then + clear_ = clear + else + clear_ = .true. + end if + + if (a%is_bld() .or. a%is_upd()) then + ! do nothing + return + else if (a%is_asb()) then + if (clear_) a%val(:) = dzero + call a%set_upd() + else + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_coo_reinit + subroutine d_coo_trim(a) use psb_realloc_mod diff --git a/base/modules/psb_d_csr_mat_mod.f03 b/base/modules/psb_d_csr_mat_mod.f03 index a67649a0..81b78e1b 100644 --- a/base/modules/psb_d_csr_mat_mod.f03 +++ b/base/modules/psb_d_csr_mat_mod.f03 @@ -29,12 +29,18 @@ module psb_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) :: mv_from => d_csr_mv_from +!!$ procedure, pass(a) :: cp_from => d_csr_cp_from + procedure, pass(a) :: csgetptn => d_csr_csgetptn procedure, pass(a) :: d_csgetrow => d_csr_csgetrow + procedure, pass(a) :: get_nz_row => d_csr_get_nz_row 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 procedure, pass(a) :: sizeof => d_csr_sizeof + procedure, pass(a) :: reinit => d_csr_reinit + end type psb_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, & @@ -44,7 +50,9 @@ module psb_d_csr_mat_mod & 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_get_size, & - & d_csr_sizeof + & d_csr_sizeof, d_csr_csgetptn, d_csr_get_nz_row, d_csr_reinit +!!$, & +!!$ & d_csr_mv_from, d_csr_mv_from interface @@ -149,6 +157,25 @@ module psb_d_csr_mat_mod end subroutine d_csr_csput_impl end interface + interface + subroutine d_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_d_csr_sparse_mat + implicit none + + class(psb_d_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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_csgetptn_impl + end interface + interface subroutine d_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) @@ -259,7 +286,7 @@ contains implicit none class(psb_d_csr_sparse_mat), intent(in) :: a integer :: res - res = a%irp(a%m+1)-1 + res = a%irp(a%get_nrows()+1)-1 end function d_csr_get_nzeros function d_csr_get_size(a) result(res) @@ -286,6 +313,26 @@ contains end function d_csr_get_size + + + function d_csr_get_nz_row(idx,a) result(res) + use psb_const_mod + implicit none + + class(psb_d_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer :: res + + res = 0 + + if ((1<=idx).and.(idx<=a%get_nrows())) then + res = a%irp(idx+1)-a%irp(idx) + end if + + end function d_csr_get_nz_row + + + !===================================== ! ! @@ -313,7 +360,8 @@ contains call psb_realloc(nz,a%ja,info) if (info == 0) call psb_realloc(nz,a%val,info) - if (info == 0) call psb_realloc(max(nz,a%m+1,a%n+1),a%irp,info) + if (info == 0) call psb_realloc(& + & max(nz,a%get_nrows()+1,a%get_ncols()+1),a%irp,info) if (info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -396,6 +444,49 @@ contains return end subroutine d_csr_csput + subroutine d_csr_csgetptn(imin,imax,a,nz,ia,ja,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(psb_d_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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_csgetptn_impl(imin,imax,a,nz,ia,ja,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_csgetptn + + 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 @@ -605,6 +696,55 @@ contains end subroutine d_csr_free + subroutine d_csr_reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_d_csr_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + + Integer :: err_act, info + character(len=20) :: name='reinit' + logical :: clear_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + + if (present(clear)) then + clear_ = clear + else + clear_ = .true. + end if + + if (a%is_bld() .or. a%is_upd()) then + ! do nothing + return + else if (a%is_asb()) then + if (clear_) a%val(:) = dzero + call a%set_upd() + else + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_csr_reinit + + subroutine d_csr_trim(a) use psb_realloc_mod use psb_error_mod @@ -1039,6 +1179,81 @@ contains end subroutine d_csr_print +!!$ +!!$ subroutine d_csr_cp_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_d_csr_sparse_mat), intent(out) :: a +!!$ class(psb_d_csr_sparse_mat), intent(inout) :: b +!!$ +!!$ +!!$ Integer :: err_act, info +!!$ character(len=20) :: name='cp_from' +!!$ logical, parameter :: debug=.false. +!!$ +!!$ call psb_erractionsave(err_act) +!!$ +!!$ info = 0 +!!$ +!!$ call a%allocate(b%get_nrows(),b%get_ncols(),b%get_nzeros()) +!!$ a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat +!!$ a%irp = b%irp +!!$ a%ja = b%ja +!!$ a%val = b%val +!!$ +!!$ if (info /= 0) goto 9999 +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 continue +!!$ call psb_erractionrestore(err_act) +!!$ +!!$ call psb_errpush(info,name) +!!$ +!!$ if (err_act /= psb_act_ret_) then +!!$ call psb_error() +!!$ end if +!!$ return +!!$ +!!$ end subroutine d_csr_cp_from +!!$ +!!$ subroutine d_csr_mv_from(a,b) +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ class(psb_d_csr_sparse_mat), intent(out) :: a +!!$ class(psb_d_csr_sparse_mat), intent(inout) :: b +!!$ +!!$ +!!$ Integer :: err_act, info +!!$ character(len=20) :: name='mv_from' +!!$ logical, parameter :: debug=.false. +!!$ +!!$ call psb_erractionsave(err_act) +!!$ info = 0 +!!$ call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat) +!!$ call move_alloc(b%irp, a%irp) +!!$ call move_alloc(b%ja, a%ja) +!!$ call move_alloc(b%val, a%val) +!!$ call b%free() +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 continue +!!$ call psb_erractionrestore(err_act) +!!$ +!!$ call psb_errpush(info,name) +!!$ +!!$ if (err_act /= psb_act_ret_) then +!!$ call psb_error() +!!$ end if +!!$ return +!!$ +!!$ end subroutine d_csr_mv_from +!!$ + !===================================== ! diff --git a/base/modules/psb_inter_desc_mod.f90 b/base/modules/psb_inter_desc_mod.f90 deleted file mode 100644 index eadb3dea..00000000 --- a/base/modules/psb_inter_desc_mod.f90 +++ /dev/null @@ -1,431 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 2.2 -!!$ (C) Copyright 2006/2007/2008 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -! -! -! package: psb_inter_descriptor_type -! Defines facilities for mapping between vectors belonging -! to different spaces. -! -module psb_inter_desc_mod - - use psb_inter_descriptor_type - use psb_descriptor_type - - interface psb_forward_map - subroutine psb_s_forward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - real(psb_spk_), intent(in) :: alpha,beta - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(out) :: y(:) - integer, intent(out) :: info - real(psb_spk_), optional :: work(:) - end subroutine psb_s_forward_map - subroutine psb_d_forward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - real(psb_dpk_), intent(in) :: alpha,beta - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(out) :: y(:) - integer, intent(out) :: info - real(psb_dpk_), optional :: work(:) - end subroutine psb_d_forward_map - subroutine psb_c_forward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - complex(psb_spk_), intent(in) :: alpha,beta - complex(psb_spk_), intent(inout) :: x(:) - complex(psb_spk_), intent(out) :: y(:) - integer, intent(out) :: info - complex(psb_spk_), optional :: work(:) - end subroutine psb_c_forward_map - subroutine psb_z_forward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - complex(psb_dpk_), intent(in) :: alpha,beta - complex(psb_dpk_), intent(inout) :: x(:) - complex(psb_dpk_), intent(out) :: y(:) - integer, intent(out) :: info - complex(psb_dpk_), optional :: work(:) - end subroutine psb_z_forward_map - end interface - - interface psb_backward_map - subroutine psb_s_backward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - real(psb_spk_), intent(in) :: alpha,beta - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(out) :: y(:) - integer, intent(out) :: info - real(psb_spk_), optional :: work(:) - end subroutine psb_s_backward_map - subroutine psb_d_backward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - real(psb_dpk_), intent(in) :: alpha,beta - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(out) :: y(:) - integer, intent(out) :: info - real(psb_dpk_), optional :: work(:) - end subroutine psb_d_backward_map - subroutine psb_c_backward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - complex(psb_spk_), intent(in) :: alpha,beta - complex(psb_spk_), intent(inout) :: x(:) - complex(psb_spk_), intent(out) :: y(:) - integer, intent(out) :: info - complex(psb_spk_), optional :: work(:) - end subroutine psb_c_backward_map - subroutine psb_z_backward_map(alpha,x,beta,y,desc,info,work) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type), intent(in) :: desc - complex(psb_dpk_), intent(in) :: alpha,beta - complex(psb_dpk_), intent(inout) :: x(:) - complex(psb_dpk_), intent(out) :: y(:) - integer, intent(out) :: info - complex(psb_dpk_), optional :: work(:) - end subroutine psb_z_backward_map - end interface - - - - interface psb_is_ok_desc - module procedure psb_is_ok_inter_desc - end interface - - interface psb_is_asb_desc - module procedure psb_is_asb_inter_desc - end interface - - interface psb_inter_desc - function psb_s_inter_desc(map_kind,desc1,desc2,map_fw,map_bk,idx_fw,idx_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_s_inter_desc - type(psb_desc_type), target :: desc1, desc2 - type(psb_sspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) - end function psb_s_inter_desc - function psb_s_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_s_inter_desc_noidx - type(psb_desc_type), target :: desc1, desc2 - type(psb_sspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind - end function psb_s_inter_desc_noidx - function psb_d_inter_desc(map_kind,desc1,desc2,map_fw,map_bk,idx_fw,idx_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_d_inter_desc - type(psb_desc_type), target :: desc1, desc2 - type(psb_dspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) - end function psb_d_inter_desc - function psb_d_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_d_inter_desc_noidx - type(psb_desc_type), target :: desc1, desc2 - type(psb_dspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind - end function psb_d_inter_desc_noidx - function psb_c_inter_desc(map_kind,desc1,desc2,map_fw,map_bk,idx_fw,idx_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_c_inter_desc - type(psb_desc_type), target :: desc1, desc2 - type(psb_cspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) - end function psb_c_inter_desc - function psb_c_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_c_inter_desc_noidx - type(psb_desc_type), target :: desc1, desc2 - type(psb_cspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind - end function psb_c_inter_desc_noidx - function psb_z_inter_desc(map_kind,desc1,desc2,map_fw,map_bk,idx_fw,idx_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_z_inter_desc - type(psb_desc_type), target :: desc1, desc2 - type(psb_zspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) - end function psb_z_inter_desc - function psb_z_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) - use psb_inter_descriptor_type - implicit none - type(psb_inter_desc_type) :: psb_z_inter_desc_noidx - type(psb_desc_type), target :: desc1, desc2 - type(psb_zspmat_type), intent(in) :: map_fw, map_bk - integer, intent(in) :: map_kind - end function psb_z_inter_desc_noidx - end interface - - interface psb_sizeof - module procedure psb_itd_sizeof,& - & psb_s_map_sizeof, psb_c_map_sizeof,& - & psb_d_map_sizeof, psb_z_map_sizeof - end interface - - interface psb_linmap - subroutine psb_s_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout) - use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & - & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof - use psb_descriptor_type, only: psb_desc_type - implicit none - real(psb_spk_), intent(in) :: alpha,beta - real(psb_spk_), intent(inout) :: x(:),y(:) - type(psb_sspmat_type), intent(in) :: a_map - type(psb_desc_type), intent(in) :: cd_xt,descin, descout - end subroutine psb_s_apply_linmap - subroutine psb_d_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout) - use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & - & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof - use psb_descriptor_type, only: psb_desc_type - implicit none - real(psb_dpk_), intent(in) :: alpha,beta - real(psb_dpk_), intent(inout) :: x(:),y(:) - type(psb_dspmat_type), intent(in) :: a_map - type(psb_desc_type), intent(in) :: cd_xt,descin, descout - end subroutine psb_d_apply_linmap - subroutine psb_c_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout) - use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & - & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof - use psb_descriptor_type, only: psb_desc_type - implicit none - complex(psb_spk_), intent(in) :: alpha,beta - complex(psb_spk_), intent(inout) :: x(:),y(:) - type(psb_cspmat_type), intent(in) :: a_map - type(psb_desc_type), intent(in) :: cd_xt,descin, descout - end subroutine psb_c_apply_linmap - subroutine psb_z_apply_linmap(alpha,x,beta,y,a_map,cd_xt,descin,descout) - use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & - & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof - use psb_descriptor_type, only: psb_desc_type - implicit none - complex(psb_dpk_), intent(in) :: alpha,beta - complex(psb_dpk_), intent(inout) :: x(:),y(:) - type(psb_zspmat_type), intent(in) :: a_map - type(psb_desc_type), intent(in) :: cd_xt,descin, descout - end subroutine psb_z_apply_linmap - end interface - -contains - - function psb_cd_get_map_kind(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - Integer :: psb_cd_get_map_kind - if (psb_is_ok_desc(desc)) then - psb_cd_get_map_kind = desc%itd_data(psb_map_kind_) - else - psb_cd_get_map_kind = -1 - end if - end function psb_cd_get_map_kind - - subroutine psb_cd_set_map_kind(map_kind,desc) - implicit none - integer, intent(in) :: map_kind - type(psb_inter_desc_type), intent(inout) :: desc - - desc%itd_data(psb_map_kind_) = map_kind - - end subroutine psb_cd_set_map_kind - - function psb_cd_get_map_data(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - Integer :: psb_cd_get_map_data - if (psb_is_ok_desc(desc)) then - psb_cd_get_map_data = desc%itd_data(psb_map_data_) - else - psb_cd_get_map_data = -1 - end if - end function psb_cd_get_map_data - - subroutine psb_cd_set_map_data(map_data,desc) - implicit none - integer, intent(in) :: map_data - type(psb_inter_desc_type), intent(inout) :: desc - - - desc%itd_data(psb_map_data_) = map_data - - end subroutine psb_cd_set_map_data - - - function psb_cd_get_fw_tmp_sz(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - Integer :: psb_cd_get_fw_tmp_sz - - psb_cd_get_fw_tmp_sz = desc%itd_data(psb_fw_tmp_sz_) - end function psb_cd_get_fw_tmp_sz - - function psb_cd_get_bk_tmp_sz(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - Integer :: psb_cd_get_bk_tmp_sz - - psb_cd_get_bk_tmp_sz = desc%itd_data(psb_bk_tmp_sz_) - end function psb_cd_get_bk_tmp_sz - - subroutine psb_cd_set_fw_tmp_sz(isz,desc) - implicit none - type(psb_inter_desc_type), intent(inout) :: desc - integer, intent(in) :: isz - - desc%itd_data(psb_fw_tmp_sz_) =isz - end subroutine psb_cd_set_fw_tmp_sz - - subroutine psb_cd_set_bk_tmp_sz(isz,desc) - implicit none - type(psb_inter_desc_type), intent(inout) :: desc - integer, intent(in) :: isz - - desc%itd_data(psb_bk_tmp_sz_) =isz - - end subroutine psb_cd_set_bk_tmp_sz - - - logical function psb_is_asb_inter_desc(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - - psb_is_asb_inter_desc = .false. - if (.not.allocated(desc%itd_data)) return - if (.not.associated(desc%desc_1)) return - if (.not.associated(desc%desc_2)) return - psb_is_asb_inter_desc = & - & psb_is_asb_desc(desc%desc_1).and.psb_is_asb_desc(desc%desc_2) - - end function psb_is_asb_inter_desc - - logical function psb_is_ok_inter_desc(desc) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - - psb_is_ok_inter_desc = .false. - if (.not.allocated(desc%itd_data)) return - select case(desc%itd_data(psb_map_data_)) - case(psb_map_integer_, psb_map_single_, psb_map_complex_,& - & psb_map_double_, psb_map_double_complex_) - ! Ok go ahead - case default - ! Since it's false so far, simply return - return - end select - if (.not.associated(desc%desc_1)) return - if (.not.associated(desc%desc_2)) return - psb_is_ok_inter_desc = & - & psb_is_ok_desc(desc%desc_1).and.psb_is_ok_desc(desc%desc_2) - - end function psb_is_ok_inter_desc - - - function psb_s_map_sizeof(map) result(val) - implicit none - type(psb_s_map_type), intent(in) :: map - integer(psb_long_int_k_) :: val - - val = 0 - val = val + psb_sizeof(map%map_fw) - val = val + psb_sizeof(map%map_bk) - - end function psb_s_map_sizeof - - function psb_d_map_sizeof(map) result(val) - implicit none - type(psb_d_map_type), intent(in) :: map - integer(psb_long_int_k_) :: val - - val = 0 - val = val + psb_sizeof(map%map_fw) - val = val + psb_sizeof(map%map_bk) - - end function psb_d_map_sizeof - - function psb_c_map_sizeof(map) result(val) - implicit none - type(psb_c_map_type), intent(in) :: map - integer(psb_long_int_k_) :: val - - val = 0 - val = val + psb_sizeof(map%map_fw) - val = val + psb_sizeof(map%map_bk) - - end function psb_c_map_sizeof - - function psb_z_map_sizeof(map) result(val) - implicit none - type(psb_z_map_type), intent(in) :: map - integer(psb_long_int_k_) :: val - - val = 0 - val = val + psb_sizeof(map%map_fw) - val = val + psb_sizeof(map%map_bk) - - end function psb_z_map_sizeof - - function psb_itd_sizeof(desc) result(val) - implicit none - type(psb_inter_desc_type), intent(in) :: desc - integer(psb_long_int_k_) :: val - - val = 0 - if (allocated(desc%itd_data)) val = val + psb_sizeof_int*size(desc%itd_data) - if (allocated(desc%exch_fw_idx)) val = val + psb_sizeof_int*size(desc%exch_fw_idx) - if (allocated(desc%exch_bk_idx)) val = val + psb_sizeof_int*size(desc%exch_bk_idx) - val = val + psb_sizeof(desc%desc_fw) - val = val + psb_sizeof(desc%desc_bk) - val = val + psb_sizeof(desc%dmap) - val = val + psb_sizeof(desc%zmap) - - end function psb_itd_sizeof - - - - -end module psb_inter_desc_mod diff --git a/base/modules/psb_inter_desc_type.f90 b/base/modules/psb_inter_desc_type.f90 deleted file mode 100644 index 2f3a3114..00000000 --- a/base/modules/psb_inter_desc_type.f90 +++ /dev/null @@ -1,89 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 2.2 -!!$ (C) Copyright 2006/2007/2008 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -! -! -! package: psb_inter_descriptor_type -! Defines facilities for mapping between vectors belonging -! to different spaces. -! -module psb_inter_descriptor_type - use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & - & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof - - use psb_descriptor_type, only: psb_desc_type - - - ! Inter-descriptor mapping data structures. - integer, parameter :: psb_map_kind_ = 1 - integer, parameter :: psb_map_data_ = 2 - integer, parameter :: psb_map_integer_ = 1 - integer, parameter :: psb_map_single_ = 2 - integer, parameter :: psb_map_double_ = 3 - integer, parameter :: psb_map_complex_ = 4 - integer, parameter :: psb_map_double_complex_ = 5 - - integer, parameter :: psb_fw_tmp_kind_ = 5 - integer, parameter :: psb_fw_tmp_sz_ = 6 - integer, parameter :: psb_bk_tmp_kind_ = 7 - integer, parameter :: psb_bk_tmp_sz_ = 8 - integer, parameter :: psb_itd_data_size_=20 - - - type psb_s_map_type - type(psb_sspmat_type) :: map_fw, map_bk - end type psb_s_map_type - - type psb_d_map_type - type(psb_dspmat_type) :: map_fw, map_bk - end type psb_d_map_type - - type psb_c_map_type - type(psb_cspmat_type) :: map_fw, map_bk - end type psb_c_map_type - - type psb_z_map_type - type(psb_zspmat_type) :: map_fw, map_bk - end type psb_z_map_type - - type psb_inter_desc_type - integer, allocatable :: itd_data(:) - type(psb_desc_type), pointer :: desc_1=>null(), desc_2=>null() - integer, allocatable :: exch_fw_idx(:), exch_bk_idx(:) - type(psb_desc_type) :: desc_fw, desc_bk - type(psb_s_map_type) :: smap - type(psb_d_map_type) :: dmap - type(psb_c_map_type) :: cmap - type(psb_z_map_type) :: zmap - end type psb_inter_desc_type - -end module psb_inter_descriptor_type - diff --git a/base/modules/psb_linmap_mod.f90 b/base/modules/psb_linmap_mod.f90 index 5389b3f6..0a2cef54 100644 --- a/base/modules/psb_linmap_mod.f90 +++ b/base/modules/psb_linmap_mod.f90 @@ -175,7 +175,7 @@ module psb_linmap_mod implicit none type(psb_dlinmap_type) :: psb_d_linmap type(psb_desc_type), target :: desc_X, desc_Y - type(psb_dspmat_type), intent(in) :: map_X2Y, map_Y2X + type(psb_d_sparse_mat), intent(in) :: map_X2Y, map_Y2X integer, intent(in) :: map_kind integer, intent(in), optional :: iaggr(:), naggr(:) end function psb_d_linmap @@ -471,6 +471,7 @@ contains end function psb_slinmap_sizeof function psb_dlinmap_sizeof(map) result(val) + use psb_d_mat_mod implicit none type(psb_dlinmap_type), intent(in) :: map integer(psb_long_int_k_) :: val @@ -544,7 +545,7 @@ contains implicit none type(psb_dlinmap_type), intent(out) :: out_map type(psb_desc_type), target :: desc_X, desc_Y - type(psb_dspmat_type), intent(in) :: map_X2Y, map_Y2X + type(psb_d_sparse_mat), intent(in) :: map_X2Y, map_Y2X integer, intent(in) :: map_kind integer, intent(in), optional :: iaggr(:), naggr(:) out_map = psb_linmap(map_kind,desc_X,desc_Y,map_X2Y,map_Y2X,iaggr,naggr) @@ -595,8 +596,9 @@ contains end subroutine psb_slinmap_transfer subroutine psb_dlinmap_transfer(mapin,mapout,info) - use psb_spmat_type + use psb_realloc_mod use psb_descriptor_type + use psb_mat_mod implicit none type(psb_dlinmap_type) :: mapin,mapout integer, intent(out) :: info diff --git a/base/modules/psb_linmap_type_mod.f90 b/base/modules/psb_linmap_type_mod.f90 index 5f2c3c69..aa9d512e 100644 --- a/base/modules/psb_linmap_type_mod.f90 +++ b/base/modules/psb_linmap_type_mod.f90 @@ -39,6 +39,7 @@ module psb_linmap_type_mod use psb_spmat_type, only : psb_sspmat_type, psb_dspmat_type, & & psb_cspmat_type, psb_zspmat_type, psb_spk_, psb_dpk_, psb_sizeof + use psb_d_mat_mod, only: psb_d_sparse_mat use psb_descriptor_type, only: psb_desc_type @@ -65,7 +66,7 @@ module psb_linmap_type_mod integer, allocatable :: itd_data(:), iaggr(:), naggr(:) type(psb_desc_type), pointer :: p_desc_X=>null(), p_desc_Y=>null() type(psb_desc_type) :: desc_X, desc_Y - type(psb_dspmat_type) :: map_X2Y, map_Y2X + type(psb_d_sparse_mat) :: map_X2Y, map_Y2X end type psb_dlinmap_type type psb_clinmap_type diff --git a/base/modules/psb_mat_mod.f03 b/base/modules/psb_mat_mod.f03 index 6e3270ee..17eefe52 100644 --- a/base/modules/psb_mat_mod.f03 +++ b/base/modules/psb_mat_mod.f03 @@ -27,6 +27,7 @@ module psb_d_mat_mod procedure, pass(a) :: get_nrows procedure, pass(a) :: get_ncols procedure, pass(a) :: get_nzeros + procedure, pass(a) :: get_nz_row procedure, pass(a) :: get_size procedure, pass(a) :: get_state procedure, pass(a) :: get_dupl @@ -48,17 +49,24 @@ module psb_d_mat_mod procedure, pass(a) :: free procedure, pass(a) :: trim procedure, pass(a) :: csput + procedure, pass(a) :: d_csgetptn procedure, pass(a) :: d_csgetrow procedure, pass(a) :: d_csgetblk - generic, public :: csget => d_csgetrow, d_csgetblk + generic, public :: csget => d_csgetptn, d_csgetrow, d_csgetblk procedure, pass(a) :: csclip procedure, pass(a) :: reall => reallocate_nz procedure, pass(a) :: get_neigh procedure, pass(a) :: d_cscnv procedure, pass(a) :: d_cscnv_ip generic, public :: cscnv => d_cscnv, d_cscnv_ip + procedure, pass(a) :: reinit procedure, pass(a) :: print => sparse_print - + procedure, pass(a) :: d_mv_from + generic, public :: mv_from => d_mv_from + procedure, pass(a) :: d_cp_from + generic, public :: cp_from => d_cp_from + + ! Computational routines procedure, pass(a) :: get_diag procedure, pass(a) :: csnmi @@ -80,17 +88,26 @@ module psb_d_mat_mod & is_unit, get_neigh, csall, csput, d_csgetrow,& & d_csgetblk, csclip, d_cscnv, d_cscnv_ip, & & reallocate_nz, free, trim, & - & sparse_print, & + & sparse_print, reinit, & & 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, get_diag + & set_unit, get_diag, get_nz_row, d_csgetptn, & + & d_mv_from, d_cp_from interface psb_sizeof module procedure d_sizeof end interface + interface psb_move_alloc + module procedure d_sparse_mat_move + end interface + + interface psb_clone + module procedure d_sparse_mat_clone + end interface + interface psb_csmm module procedure d_csmm, d_csmv end interface @@ -389,6 +406,22 @@ contains end function get_size + function get_nz_row(idx,a) result(res) + use psb_error_mod + implicit none + integer, intent(in) :: idx + class(psb_d_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act + + res = 0 + + if (allocated(a%a)) res = a%a%get_nz_row(idx) + + end function get_nz_row + + !===================================== ! @@ -823,8 +856,6 @@ contains !===================================== - - subroutine sparse_print(iout,a,iv,eirs,eics,head,ivr,ivc) use psb_error_mod implicit none @@ -989,7 +1020,7 @@ contains endif call a%a%free() - + deallocate(a%a) return 9999 continue @@ -1071,6 +1102,54 @@ contains end subroutine csput + subroutine d_csgetptn(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod + implicit none + + class(psb_d_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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,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_csgetptn + subroutine d_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) ! Output is always in COO format @@ -1407,7 +1486,132 @@ contains end subroutine d_cscnv_ip + subroutine d_mv_from(a,b) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_d_sparse_mat), intent(out) :: a + class(psb_d_base_sparse_mat), intent(inout) :: b + integer :: info + + allocate(a%a,source=b, stat=info) + call a%a%mv_from_fmt(b,info) + + return + end subroutine d_mv_from + + subroutine d_cp_from(a,b) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_d_sparse_mat), intent(out) :: a + class(psb_d_base_sparse_mat), intent(inout), allocatable :: b + Integer :: err_act, info + character(len=20) :: name='clone' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + allocate(a%a,source=b,stat=info) + if (info /= 0) info = 4000 + if (info == 0) call a%a%cp_from_fmt(b, 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 + end subroutine d_cp_from + + subroutine d_sparse_mat_move(a,b,info) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_d_sparse_mat), intent(inout) :: a + class(psb_d_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='move_alloc' + logical, parameter :: debug=.false. + + info = 0 + call move_alloc(a%a,b%a) + + return + end subroutine d_sparse_mat_move + + subroutine d_sparse_mat_clone(a,b,info) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_d_sparse_mat), intent(in) :: a + class(psb_d_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='clone' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + allocate(b%a,source=a%a,stat=info) + if (info /= 0) info = 4000 + if (info == 0) call b%a%cp_from_fmt(a%a, 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 + + end subroutine d_sparse_mat_clone + + + subroutine reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_d_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + Integer :: err_act, info + character(len=20) :: name='reinit' + + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%reinit(clear) + + 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 reinit !===================================== @@ -1554,7 +1758,7 @@ contains call psb_errpush(info,name) goto 9999 endif - + call a%a%cssm(alpha,x,beta,y,info,trans,side,d) if (info /= 0) goto 9999 @@ -1720,3 +1924,6 @@ contains end module psb_d_mat_mod +module psb_mat_mod + use psb_d_mat_mod +end module psb_mat_mod diff --git a/base/modules/psb_spmat_type.f03 b/base/modules/psb_spmat_type.f03 index da6f3d9a..ecf99412 100644 --- a/base/modules/psb_spmat_type.f03 +++ b/base/modules/psb_spmat_type.f03 @@ -115,13 +115,13 @@ module psb_spmat_type type, extends(psb_base_spmat_type) :: psb_dspmat_type real(psb_dpk_), allocatable :: aspk(:) - contains - procedure, pass(a) :: psb_dcsmm - procedure, pass(a) :: psb_dcsmv - generic, public :: csmm => psb_dcsmm, psb_dcsmv - procedure, pass(t) :: psb_dcssm - procedure, pass(t) :: psb_dcssv - generic, public :: cssm => psb_dcssm, psb_dcssv +!!$ contains +!!$ procedure, pass(a) :: psb_dcsmm +!!$ procedure, pass(a) :: psb_dcsmv +!!$ generic, public :: csmm => psb_dcsmm, psb_dcsmv +!!$ procedure, pass(t) :: psb_dcssm +!!$ procedure, pass(t) :: psb_dcssv +!!$ generic, public :: cssm => psb_dcssm, psb_dcssv end type psb_dspmat_type type, extends(psb_base_spmat_type) :: psb_zspmat_type diff --git a/base/modules/psb_tools_mod.f90 b/base/modules/psb_tools_mod.f90 index 93adca3c..d2b65680 100644 --- a/base/modules/psb_tools_mod.f90 +++ b/base/modules/psb_tools_mod.f90 @@ -33,7 +33,6 @@ Module psb_tools_mod use psb_const_mod use psb_spmat_type - interface psb_cd_set_bld subroutine psb_cd_set_bld(desc,info) use psb_descriptor_type @@ -232,9 +231,9 @@ Module psb_tools_mod Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& & rowscale,colscale,outfmt,data) use psb_descriptor_type - use psb_spmat_type - Type(psb_dspmat_type),Intent(in) :: a - Type(psb_dspmat_type),Intent(inout) :: blk + use psb_mat_mod + Type(psb_d_sparse_mat),Intent(in) :: a + Type(psb_d_sparse_mat),Intent(inout) :: blk Type(psb_desc_type),Intent(in),target :: desc_a integer, intent(out) :: info logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale @@ -488,9 +487,9 @@ Module psb_tools_mod end Subroutine psb_scdbldext Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info,extype) use psb_descriptor_type - Use psb_spmat_type + Use psb_mat_mod integer, intent(in) :: novr - Type(psb_dspmat_type), Intent(in) :: a + Type(psb_d_sparse_mat), Intent(in) :: a Type(psb_desc_type), Intent(in), target :: desc_a Type(psb_desc_type), Intent(out) :: desc_ov integer, intent(out) :: info @@ -686,10 +685,10 @@ Module psb_tools_mod end subroutine psb_dspins subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) use psb_descriptor_type - use psb_spmat_type + use psb_d_mat_mod + type(psb_d_sparse_mat), intent(inout) :: a type(psb_desc_type), intent(in) :: desc_ar type(psb_desc_type), intent(inout) :: desc_ac - type(psb_dspmat_type), intent(inout) :: a integer, intent(in) :: nz,ia(:),ja(:) real(psb_dpk_), intent(in) :: val(:) integer, intent(out) :: info @@ -748,9 +747,9 @@ Module psb_tools_mod end subroutine psb_ssprn subroutine psb_dsprn(a, desc_a,info,clear) use psb_descriptor_type - use psb_spmat_type + use psb_mat_mod type(psb_desc_type), intent(in) :: desc_a - type(psb_dspmat_type), intent(inout) :: a + type(psb_d_sparse_mat), intent(inout) :: a integer, intent(out) :: info logical, intent(in), optional :: clear end subroutine psb_dsprn @@ -1159,13 +1158,13 @@ contains subroutine psb_dlinmap_init(a_map,cd_xt,descin,descout) - use psb_spmat_type use psb_descriptor_type use psb_serial_mod use psb_penv_mod use psb_error_mod + use psb_mat_mod implicit none - type(psb_dspmat_type), intent(out) :: a_map + type(psb_d_sparse_mat), intent(out) :: a_map type(psb_desc_type), intent(out) :: cd_xt type(psb_desc_type), intent(in) :: descin, descout @@ -1185,18 +1184,18 @@ contains ncol_in = psb_cd_get_local_cols(cd_xt) nrow_out = psb_cd_get_local_rows(descout) - call psb_sp_all(nrow_out,ncol_in,a_map,info) + call a_map%csall(nrow_out,ncol_in,info) end subroutine psb_dlinmap_init subroutine psb_dlinmap_ins(nz,ir,ic,val,a_map,cd_xt,descin,descout) - use psb_spmat_type + use psb_mat_mod use psb_descriptor_type implicit none integer, intent(in) :: nz integer, intent(in) :: ir(:),ic(:) real(kind(1.d0)), intent(in) :: val(:) - type(psb_dspmat_type), intent(inout) :: a_map + type(psb_d_sparse_mat), intent(inout) :: a_map type(psb_desc_type), intent(inout) :: cd_xt type(psb_desc_type), intent(in) :: descin, descout integer :: info @@ -1205,11 +1204,11 @@ contains end subroutine psb_dlinmap_ins subroutine psb_dlinmap_asb(a_map,cd_xt,descin,descout,afmt) - use psb_spmat_type + use psb_mat_mod use psb_descriptor_type use psb_serial_mod implicit none - type(psb_dspmat_type), intent(inout) :: a_map + type(psb_d_sparse_mat), intent(inout) :: a_map type(psb_desc_type), intent(inout) :: cd_xt type(psb_desc_type), intent(in) :: descin, descout character(len=*), optional, intent(in) :: afmt @@ -1220,8 +1219,8 @@ contains ictxt = psb_cd_get_context(descin) call psb_cdasb(cd_xt,info) - a_map%k = psb_cd_get_local_cols(cd_xt) - call psb_spcnv(a_map,info,afmt=afmt) + call a_map%set_ncols(psb_cd_get_local_cols(cd_xt)) + call a_map%cscnv(info,type=afmt) end subroutine psb_dlinmap_asb diff --git a/base/psblas/psb_dnrmi.f90 b/base/psblas/psb_dnrmi.f90 index 41e6c758..ea75ae7c 100644 --- a/base/psblas/psb_dnrmi.f90 +++ b/base/psblas/psb_dnrmi.f90 @@ -43,17 +43,16 @@ ! function psb_dnrmi(a,desc_a,info) use psb_descriptor_type - use psb_serial_mod use psb_check_mod use psb_error_mod use psb_penv_mod use psb_d_mat_mod implicit none - type(psb_d_sparse_mat), intent(in) :: a - integer, intent(out) :: info - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_) :: psb_dnrmi + type(psb_d_sparse_mat), intent(in) :: a + integer, intent(out) :: info + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_) :: psb_dnrmi ! locals integer :: ictxt, np, me,& diff --git a/base/psblas/psb_dspmm.f90 b/base/psblas/psb_dspmm.f90 index cfe50acf..fe5a82cf 100644 --- a/base/psblas/psb_dspmm.f90 +++ b/base/psblas/psb_dspmm.f90 @@ -65,8 +65,6 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& & trans, k, jx, jy, work, doswap) - use psb_spmat_type - use psb_serial_mod use psb_descriptor_type use psb_comm_mod use psi_mod @@ -426,8 +424,6 @@ end subroutine psb_dspmm subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& & trans, work, doswap) - use psb_spmat_type - use psb_serial_mod use psb_descriptor_type use psb_comm_mod use psb_const_mod @@ -443,7 +439,6 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& real(psb_dpk_), intent(inout), target :: x(:) real(psb_dpk_), intent(inout), target :: y(:) type(psb_d_sparse_mat), intent(in) :: a -!!$ type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(psb_dpk_), optional, target :: work(:) diff --git a/base/psblas/psb_dspsm.f90 b/base/psblas/psb_dspsm.f90 index 262249df..1fe76546 100644 --- a/base/psblas/psb_dspsm.f90 +++ b/base/psblas/psb_dspsm.f90 @@ -77,8 +77,6 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& & trans, side, choice, diag, k, jx, jy, work) - use psb_spmat_type - use psb_serial_mod use psb_descriptor_type use psb_comm_mod use psi_mod @@ -364,8 +362,6 @@ end subroutine psb_dspsm ! subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& & trans, side, choice, diag, work) - use psb_spmat_type - use psb_serial_mod use psb_descriptor_type use psb_comm_mod use psi_mod diff --git a/base/serial/Makefile b/base/serial/Makefile index accd7084..8913e629 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -1,20 +1,16 @@ include ../../Make.inc -FOBJS = psb_cest.o psb_dcoins.o psb_dcsmm.o psb_dcsmv.o \ - psb_dcsnmi.o psb_dcsprt.o psb_dcsrws.o psb_dcssm.o psb_dcssv.o \ - psb_dfixcoo.o psb_dipcoo2csr.o psb_dipcsr2coo.o psb_dneigh.o \ - psb_dnumbmm.o psb_drwextd.o psb_dspgtdiag.o psb_dspgtblk.o \ - psb_dspscal.o psb_dsymbmm.o psb_dtransp.o psb_dspclip.o psb_dspcnv.o\ - psb_regen_mod.o psb_dipcoo2csc.o psb_dspgetrow.o psb_lsame.o psb_zspgetrow.o\ +FOBJS = psb_cest.o \ + psb_regen_mod.o psb_lsame.o psb_zspgetrow.o\ psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtblk.o\ psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\ psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \ psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o psb_ztransc.o\ psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspscal.o psb_zspclip.o\ psb_getifield.o psb_setifield.o psb_update_mod.o psb_getrow_mod.o\ - psb_dgelp.o psb_zgelp.o\ - psb_dspshift.o psb_dspsetbld.o psb_zspshift.o psb_zspsetbld.o\ + psb_zgelp.o\ + psb_zspshift.o psb_zspsetbld.o\ psb_scsprt.o psb_sspcnv.o psb_scoins.o psb_scsmm.o psb_scsmv.o \ psb_scssm.o psb_scssv.o psb_sneigh.o psb_sspgtblk.o psb_sspgetrow.o \ psb_sfixcoo.o psb_sipcsr2coo.o psb_sipcoo2csr.o psb_sipcoo2csc.o \ @@ -27,6 +23,32 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsmm.o psb_dcsmv.o \ psb_ccssm.o psb_ccssv.o psb_ccsmm.o psb_ccsmv.o psb_ctransp.o psb_ctransc.o\ psb_cspclip.o psb_crwextd.o psb_cspscal.o\ psb_cnumbmm.o psb_csymbmm.o psb_cneigh.o +#FOBJS = psb_cest.o psb_dcoins.o psb_dcsmm.o psb_dcsmv.o \ +# psb_dcsnmi.o psb_dcsprt.o psb_dcsrws.o psb_dcssm.o psb_dcssv.o \ +# psb_dfixcoo.o psb_dipcoo2csr.o psb_dipcsr2coo.o psb_dneigh.o \ +# psb_dnumbmm.o psb_drwextd.o psb_dspgtdiag.o psb_dspgtblk.o \ +# psb_dspscal.o psb_dsymbmm.o psb_dtransp.o psb_dspclip.o psb_dspcnv.o\ +# psb_regen_mod.o psb_dipcoo2csc.o psb_dspgetrow.o psb_lsame.o psb_zspgetrow.o\ +# psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtblk.o\ +# psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\ +# psb_zfixcoo.o psb_zipcoo2csr.o psb_zipcsr2coo.o psb_zipcoo2csc.o \ +# psb_zcoins.o psb_zcsprt.o psb_zneigh.o psb_ztransp.o psb_ztransc.o\ +# psb_zrwextd.o psb_zsymbmm.o psb_znumbmm.o psb_zspscal.o psb_zspclip.o\ +# psb_getifield.o psb_setifield.o psb_update_mod.o psb_getrow_mod.o\ +# psb_dgelp.o psb_zgelp.o\ +# psb_dspshift.o psb_dspsetbld.o psb_zspshift.o psb_zspsetbld.o\ +# psb_scsprt.o psb_sspcnv.o psb_scoins.o psb_scsmm.o psb_scsmv.o \ +# psb_scssm.o psb_scssv.o psb_sneigh.o psb_sspgtblk.o psb_sspgetrow.o \ +# psb_sfixcoo.o psb_sipcsr2coo.o psb_sipcoo2csr.o psb_sipcoo2csc.o \ +# psb_sgelp.o psb_sspgtdiag.o psb_scsnmi.o psb_stransp.o \ +# psb_sspclip.o psb_srwextd.o psb_sspscal.o\ +# psb_snumbmm.o psb_ssymbmm.o\ +# psb_ccsprt.o psb_cspcnv.o psb_ccoins.o psb_ccsnmi.o\ +# psb_cfixcoo.o psb_cipcsr2coo.o psb_cipcoo2csr.o psb_cipcoo2csc.o \ +# psb_cgelp.o psb_cspgtdiag.o psb_cspgtblk.o psb_cspgetrow.o\ +# psb_ccssm.o psb_ccssv.o psb_ccsmm.o psb_ccsmv.o psb_ctransp.o psb_ctransc.o\ +# psb_cspclip.o psb_crwextd.o psb_cspscal.o\ +# psb_cnumbmm.o psb_csymbmm.o psb_cneigh.o # psb_dcsrp.o psb_zcsrp.o\ LIBDIR=.. diff --git a/base/serial/coo/Makefile b/base/serial/coo/Makefile index 84bdcb23..ae35f209 100644 --- a/base/serial/coo/Makefile +++ b/base/serial/coo/Makefile @@ -4,10 +4,10 @@ include ../../../Make.inc # The object files # FOBJS = scoonrmi.o scoomm.o scoomv.o scoosm.o scoosv.o scoorws.o\ - dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\ ccoonrmi.o ccoomm.o ccoomv.o ccoosm.o ccoosv.o ccoorws.o\ zcoomm.o zcoomv.o zcoonrmi.o zcoorws.o zcoosm.o zcoosv.o +# dcoonrmi.o dcoomm.o dcoomv.o dcoosm.o dcoosv.o dcoorws.o\ OBJS=$(FOBJS) diff --git a/base/serial/csr/Makefile b/base/serial/csr/Makefile index 6349bd3f..51daae9f 100644 --- a/base/serial/csr/Makefile +++ b/base/serial/csr/Makefile @@ -4,13 +4,14 @@ include ../../../Make.inc # The object files # -FOBJS = dcsrck.o dcsrmm.o dcsrsm.o dcsrmv.o dcsrsv.o dcrnrmi.o \ - dcsrmv4.o dcsrmv2.o dcsrmv3.o dcsrrws.o\ - scsrmm.o scsrmv.o scsrmv2.o scsrmv3.o scsrmv4.o scsrsm.o scsrsv.o\ +FOBJS = scsrmm.o scsrmv.o scsrmv2.o scsrmv3.o scsrmv4.o scsrsm.o scsrsv.o\ scrnrmi.o \ ccrnrmi.o ccsrmm.o ccsrrws.o ccsrsm.o csrmv.o csrsv.o ccsrck.o\ zcrnrmi.o zcsrmm.o zcsrrws.o zcsrsm.o zsrmv.o zsrsv.o zcsrck.o +#dcsrck.o dcsrmm.o dcsrsm.o dcsrmv.o dcsrsv.o dcrnrmi.o \ +# dcsrmv4.o dcsrmv2.o dcsrmv3.o dcsrrws.o\ + OBJS=$(FOBJS) # diff --git a/base/serial/dp/Makefile b/base/serial/dp/Makefile index 7594c686..b175cb51 100644 --- a/base/serial/dp/Makefile +++ b/base/serial/dp/Makefile @@ -3,16 +3,21 @@ include ../../../Make.inc # The object files # -XOBJS = scrjd.o dcrjd.o ccrjd.o zcrjd.o -FOBJS = dcrcr.o dgblock.o partition.o \ - dgindex.o djadrp.o djadrp1.o dcsrrp.o dcsrp1.o check_dim.o \ - Max_nnzero.o dcoco.o dcocr.o dcrco.o djdcox.o djdco.o dvtfg.o dgind_tri.o \ +XOBJS = scrjd.o ccrjd.o zcrjd.o +FOBJS = partition.o dgblock.o dvtfg.o \ + check_dim.o \ + Max_nnzero.o \ gen_block.o\ scrco.o scrcr.o scocr.o scoco.o sgindex.o sgind_tri.o\ ccoco.o ccocr.o ccrco.o ccrcr.o cgindex.o cgind_tri.o\ zcoco.o zcocr.o zcrco.o zcrcr.o zgindex.o zgind_tri.o\ $(XOBJS) +#dcoco.o dcocr.o dcrco.o djdcox.o djdco.o dgind_tri.o \ +#dcrcr.o +#dgindex.o djadrp.o djadrp1.o dcsrrp.o dcsrp1.o +#dcrjd.o + # # dgind_tri.o # diff --git a/base/serial/f03/psbn_d_coo_impl.f03 b/base/serial/f03/psbn_d_coo_impl.f03 index d696675c..e9bb6b48 100644 --- a/base/serial/f03/psbn_d_coo_impl.f03 +++ b/base/serial/f03/psbn_d_coo_impl.f03 @@ -902,6 +902,277 @@ end function d_coo_csnmi_impl +subroutine d_coo_csgetptn_impl(imin,imax,a,nz,ia,ja,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 psb_d_base_mat_mod, psb_protect_name => d_coo_csgetptn_impl + implicit none + + class(psb_d_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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) 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 + 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 + 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) 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) return + end if + 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) return + + end if + 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_getptn + +end subroutine d_coo_csgetptn_impl + + 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 @@ -1352,7 +1623,7 @@ contains integer, intent(out) :: info integer, intent(in), optional :: gtl(:) integer :: i,ir,ic, ilr, ilc, ip, & - & i1,i2,nc,nnz,dupl,ng + & i1,i2,nc,nnz,dupl,ng, nr integer :: debug_level, debug_unit character(len=20) :: name='d_coo_srch_upd' @@ -1370,6 +1641,8 @@ contains ilr = -1 ilc = -1 nnz = a%get_nzeros() + nr = a%get_nrows() + nc = a%get_ncols() if (present(gtl)) then @@ -1384,7 +1657,7 @@ contains ic = ja(i) if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then ir = gtl(ir) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then ic = gtl(ic) if (ir /= ilr) then i1 = psb_ibsrch(ir,nnz,a%ia) @@ -1427,7 +1700,7 @@ contains if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then ir = gtl(ir) ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then if (ir /= ilr) then i1 = psb_ibsrch(ir,nnz,a%ia) @@ -1479,7 +1752,7 @@ contains do i=1, nz ir = ia(i) ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then if (ir /= ilr) then i1 = psb_ibsrch(ir,nnz,a%ia) @@ -1515,7 +1788,7 @@ contains do i=1, nz ir = ia(i) ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then if (ir /= ilr) then i1 = psb_ibsrch(ir,nnz,a%ia) @@ -1576,15 +1849,9 @@ subroutine d_cp_coo_to_coo_impl(a,b,info) call psb_erractionsave(err_act) info = 0 - call b%set_nzeros(a%get_nzeros()) - call b%set_nrows(a%get_nrows()) - call b%set_ncols(a%get_ncols()) - call b%set_dupl(a%get_dupl()) - call b%set_state(a%get_state()) - call b%set_triangle(a%is_triangle()) - call b%set_upper(a%is_upper()) - call b%set_unit(a%is_unit()) + b%psb_d_base_sparse_mat = a%psb_d_base_sparse_mat + call b%set_nzeros(a%get_nzeros()) call b%reallocate(a%get_nzeros()) b%ia(:) = a%ia(:) @@ -1615,7 +1882,7 @@ subroutine d_cp_coo_from_coo_impl(a,b,info) use psb_realloc_mod use psb_d_base_mat_mod, psb_protect_name => d_cp_coo_from_coo_impl implicit none - class(psb_d_coo_sparse_mat), intent(inout) :: a + class(psb_d_coo_sparse_mat), intent(out) :: a class(psb_d_coo_sparse_mat), intent(in) :: b integer, intent(out) :: info @@ -1627,15 +1894,8 @@ subroutine d_cp_coo_from_coo_impl(a,b,info) call psb_erractionsave(err_act) info = 0 + a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat call a%set_nzeros(b%get_nzeros()) - call a%set_nrows(b%get_nrows()) - call a%set_ncols(b%get_ncols()) - call a%set_dupl(b%get_dupl()) - call a%set_state(b%get_state()) - call a%set_triangle(b%is_triangle()) - call a%set_upper(b%is_upper()) - call a%set_unit(b%is_unit()) - call a%reallocate(b%get_nzeros()) a%ia(:) = b%ia(:) @@ -2051,15 +2311,8 @@ subroutine d_mv_coo_to_coo_impl(a,b,info) call psb_erractionsave(err_act) info = 0 + b%psb_d_base_sparse_mat = a%psb_d_base_sparse_mat call b%set_nzeros(a%get_nzeros()) - call b%set_nrows(a%get_nrows()) - call b%set_ncols(a%get_ncols()) - call b%set_dupl(a%get_dupl()) - call b%set_state(a%get_state()) - call b%set_triangle(a%is_triangle()) - call b%set_upper(a%is_upper()) - call b%set_unit(a%is_unit()) - call b%reallocate(a%get_nzeros()) call move_alloc(a%ia, b%ia) @@ -2103,15 +2356,8 @@ subroutine d_mv_coo_from_coo_impl(a,b,info) call psb_erractionsave(err_act) info = 0 + a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat call a%set_nzeros(b%get_nzeros()) - call a%set_nrows(b%get_nrows()) - call a%set_ncols(b%get_ncols()) - call a%set_dupl(b%get_dupl()) - call a%set_state(b%get_state()) - call a%set_triangle(b%is_triangle()) - call a%set_upper(b%is_upper()) - call a%set_unit(b%is_unit()) - call a%reallocate(b%get_nzeros()) call move_alloc(b%ia , a%ia ) diff --git a/base/serial/f03/psbn_d_csr_impl.f03 b/base/serial/f03/psbn_d_csr_impl.f03 index 840e4cf7..89bdaed7 100644 --- a/base/serial/f03/psbn_d_csr_impl.f03 +++ b/base/serial/f03/psbn_d_csr_impl.f03 @@ -1046,6 +1046,178 @@ end function d_csr_csnmi_impl !===================================== +subroutine d_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,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 psb_d_base_mat_mod + use psb_d_csr_mat_mod, psb_protect_name => d_csr_csgetptn_impl + implicit none + + class(psb_d_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + 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=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then ir = gtl(ir) ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then i1 = a%irp(ir) i2 = a%irp(ir+1) nc=i2-i1 @@ -1370,7 +1544,7 @@ contains if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then ir = gtl(ir) ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then i1 = a%irp(ir) i2 = a%irp(ir+1) nc = i2-i1 @@ -1414,7 +1588,7 @@ contains ir = ia(i) ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then i1 = a%irp(ir) i2 = a%irp(ir+1) @@ -1447,7 +1621,7 @@ contains do i=1, nz ir = ia(i) ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then + if ((ir > 0).and.(ir <= nr)) then i1 = a%irp(ir) i2 = a%irp(ir+1) nc = i2-i1 @@ -1534,6 +1708,7 @@ subroutine d_cp_csr_to_coo_impl(a,b,info) nza = a%get_nzeros() call b%allocate(nr,nc,nza) + b%psb_d_base_sparse_mat = a%psb_d_base_sparse_mat do i=1, nr do j=a%irp(i),a%irp(i+1)-1 @@ -1542,15 +1717,7 @@ subroutine d_cp_csr_to_coo_impl(a,b,info) b%val(j) = a%val(j) end do end do - call b%set_nzeros(a%get_nzeros()) - call b%set_nrows(a%get_nrows()) - call b%set_ncols(a%get_ncols()) - call b%set_dupl(a%get_dupl()) - call b%set_state(a%get_state()) - call b%set_triangle(a%is_triangle()) - call b%set_upper(a%is_upper()) - call b%set_unit(a%is_unit()) call b%fix(info) @@ -1582,15 +1749,9 @@ subroutine d_mv_csr_to_coo_impl(a,b,info) nc = a%get_ncols() nza = a%get_nzeros() - call b%set_nzeros(a%get_nzeros()) - call b%set_nrows(a%get_nrows()) - call b%set_ncols(a%get_ncols()) - call b%set_dupl(a%get_dupl()) - call b%set_state(a%get_state()) - call b%set_triangle(a%is_triangle()) - call b%set_upper(a%is_upper()) - call b%set_unit(a%is_unit()) + b%psb_d_base_sparse_mat = a%psb_d_base_sparse_mat + call b%set_nzeros(a%get_nzeros()) call move_alloc(a%ja,b%ja) call move_alloc(a%val,b%val) call psb_realloc(nza,b%ia,info) @@ -1635,14 +1796,8 @@ subroutine d_mv_csr_from_coo_impl(a,b,info) nr = b%get_nrows() nc = b%get_ncols() nza = b%get_nzeros() - - call a%set_nrows(b%get_nrows()) - call a%set_ncols(b%get_ncols()) - call a%set_dupl(b%get_dupl()) - call a%set_state(b%get_state()) - call a%set_triangle(b%is_triangle()) - call a%set_upper(b%is_upper()) - call a%set_unit(b%is_unit()) + + a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat ! Dirty trick: call move_alloc to have the new data allocated just once. call move_alloc(b%ia,itemp) call move_alloc(b%ja,a%ja) @@ -1725,11 +1880,16 @@ subroutine d_mv_csr_to_fmt_impl(a,b,info) info = 0 select type (b) - class is (psb_d_coo_sparse_mat) + type is (psb_d_coo_sparse_mat) call a%mv_to_coo(b,info) ! Need to fix trivial copies! -!!$ class is (psb_d_csr_sparse_mat) -!!$ call a%mv_to_coo(b,info) + type is (psb_d_csr_sparse_mat) + b%psb_d_base_sparse_mat = a%psb_d_base_sparse_mat + call move_alloc(a%irp, b%irp) + call move_alloc(a%ja, b%ja) + call move_alloc(a%val, b%val) + call a%free() + class default call tmp%mv_from_fmt(a,info) if (info == 0) call b%mv_from_coo(tmp,info) @@ -1761,8 +1921,12 @@ subroutine d_cp_csr_to_fmt_impl(a,b,info) select type (b) - class is (psb_d_coo_sparse_mat) + type is (psb_d_coo_sparse_mat) call a%cp_to_coo(b,info) + + type is (psb_d_csr_sparse_mat) + b = a + class default call tmp%cp_from_fmt(a,info) if (info == 0) call b%mv_from_coo(tmp,info) @@ -1793,8 +1957,16 @@ subroutine d_mv_csr_from_fmt_impl(a,b,info) info = 0 select type (b) - class is (psb_d_coo_sparse_mat) + type is (psb_d_coo_sparse_mat) call a%mv_from_coo(b,info) + + type is (psb_d_csr_sparse_mat) + a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat + call move_alloc(b%irp, a%irp) + call move_alloc(b%ja, a%ja) + call move_alloc(b%val, a%val) + call b%free() + class default call tmp%mv_from_fmt(b,info) if (info == 0) call a%mv_from_coo(tmp,info) @@ -1818,7 +1990,7 @@ subroutine d_cp_csr_from_fmt_impl(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer :: nz, nr, i,j,irw, idl,err_act, nc Integer, Parameter :: maxtry=8 integer :: debug_level, debug_unit character(len=20) :: name @@ -1826,8 +1998,15 @@ subroutine d_cp_csr_from_fmt_impl(a,b,info) info = 0 select type (b) - class is (psb_d_coo_sparse_mat) + type is (psb_d_coo_sparse_mat) call a%cp_from_coo(b,info) + + type is (psb_d_csr_sparse_mat) + a%psb_d_base_sparse_mat = b%psb_d_base_sparse_mat + a%irp = b%irp + a%ja = b%ja + a%val = b%val + class default call tmp%cp_from_fmt(b,info) if (info == 0) call a%mv_from_coo(tmp,info) diff --git a/base/serial/f77/Makefile b/base/serial/f77/Makefile index 2fc6fbcb..e4c47c89 100644 --- a/base/serial/f77/Makefile +++ b/base/serial/f77/Makefile @@ -3,16 +3,18 @@ include ../../../Make.inc # # The object files # -FOBJS = daxpby.o dcsmm.o dcsnmi.o dcsrp.o dcssm.o \ - dgelp.o dlpupd.o dswmm.o \ - dswsm.o smmp.o dcsrws.o \ - saxpby.o slpupd.o scsmm.o sswmm.o scsnmi.o scsrws.o\ +FOBJS = daxpby.o saxpby.o slpupd.o scsmm.o sswmm.o scsnmi.o scsrws.o\ sswsm.o scssm.o sgelp.o\ caxpby.o clpupd.o ccsmm.o cswmm.o ccsnmi.o ccsrws.o\ cswsm.o ccssm.o cgelp.o\ zcsnmi.o zaxpby.o zcsmm.o zcssm.o zswmm.o zswsm.o\ zcsrws.o zgelp.o zlpupd.o +#dcsmm.o dcsnmi.o dcsrp.o dcssm.o \ +# dgelp.o dlpupd.o dswmm.o \ +# dswsm.o smmp.o dcsrws.o \ + + OBJS=$(FOBJS) # diff --git a/base/serial/jad/Makefile b/base/serial/jad/Makefile index d4ad4396..5e664ab6 100644 --- a/base/serial/jad/Makefile +++ b/base/serial/jad/Makefile @@ -3,11 +3,13 @@ include ../../../Make.inc # The object files # -FOBJS = djadmm.o djadmv.o djadsm.o djadsv.o djdnrmi.o djadnr.o\ - djadmv2.o djadmv3.o djadmv4.o djadrws.o djdrws.o \ - sjadmm.o sjadmv.o sjadsm.o sjadsv.o sjdnrmi.o sjadnr.o\ +FOBJS = sjadmm.o sjadmv.o sjadsm.o sjadsv.o sjdnrmi.o sjadnr.o\ sjadmv2.o sjadmv3.o sjadmv4.o sjadrws.o sjdrws.o +#djadmm.o djadmv.o djadsm.o djadsv.o djdnrmi.o djadnr.o\ +# djadmv2.o djadmv3.o djadmv4.o djadrws.o djdrws.o \ + + OBJS=$(FOBJS) # diff --git a/base/serial/psb_getrow_mod.f90 b/base/serial/psb_getrow_mod.f90 index ad771f54..ea1e343b 100644 --- a/base/serial/psb_getrow_mod.f90 +++ b/base/serial/psb_getrow_mod.f90 @@ -32,13 +32,13 @@ module psb_getrow_mod interface csr_getrow - module procedure csr_sspgtrow, csr_dspgtrow, csr_cspgtrow, csr_zspgtrow + module procedure csr_sspgtrow, csr_cspgtrow, csr_zspgtrow end interface interface coo_getrow - module procedure coo_sspgtrow, coo_dspgtrow, coo_cspgtrow, coo_zspgtrow + module procedure coo_sspgtrow, coo_cspgtrow, coo_zspgtrow end interface interface jad_getrow - module procedure jad_sspgtrow, jad_dspgtrow, jad_cspgtrow, jad_zspgtrow + module procedure jad_sspgtrow, jad_cspgtrow, jad_zspgtrow end interface contains @@ -468,429 +468,429 @@ contains end subroutine jad_sspgtrow - subroutine csr_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) - - use psb_sort_mod - use psb_spmat_type - use psb_const_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - integer :: irw - integer, intent(out) :: nz - integer, allocatable, intent(inout) :: ia(:), ja(:) - real(psb_dpk_), allocatable, intent(inout) :: val(:) - integer :: nzin - logical, intent(in) :: append - integer :: lrw,info - integer, optional :: iren(:) - - integer :: idx,i,j, k, nr, row_idx, nzin_ - integer, allocatable :: indices(:) - - if (append) then - nzin_ = nzin - else - nzin_ = 0 - endif - - if (a%pl(1) /= 0) then - - nr = lrw - irw + 1 - allocate(indices(nr)) - nz = 0 - do i=1,nr - indices(i)=a%pl(irw+i-1) - nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i)) - end do - - call psb_ensure_size(nzin_+nz,ia,info) - if (info==0) call psb_ensure_size(nzin_+nz,ja,info) - if (info==0) call psb_ensure_size(nzin_+nz,val,info) - if (info /= 0) return - - k=0 - if(present(iren)) then - do i=1,nr - row_idx=indices(i) - do j=a%ia2(row_idx),a%ia2(row_idx+1)-1 - k = k + 1 - val(nzin_+k) = a%aspk(j) - ia(nzin_+k) = iren(row_idx) - ja(nzin_+k) = iren(a%ia1(j)) - end do - end do - else - do i=1,nr - row_idx=indices(i) - do j=a%ia2(row_idx),a%ia2(row_idx+1)-1 - k = k + 1 - val(nzin_+k) = a%aspk(j) - ia(nzin_+k) = row_idx - ja(nzin_+k) = a%ia1(j) - end do - end do - end if - - else - - idx = irw - - if (idx<0) then - write(0,*) ' spgtrow Error : idx no good ',idx - info = 2 - return - end if - nr = lrw - irw + 1 - nz = a%ia2(idx+nr) - a%ia2(idx) - - call psb_ensure_size(nzin_+nz,ia,info) - if (info==0) call psb_ensure_size(nzin_+nz,ja,info) - if (info==0) call psb_ensure_size(nzin_+nz,val,info) - if (info /= 0) return - - - if (present(iren)) then - k=0 - do i=irw,lrw - do j=a%ia2(i),a%ia2(i+1)-1 - k = k + 1 - val(nzin_+k) = a%aspk(j) - ia(nzin_+k) = iren(i) - ja(nzin_+k) = iren(a%ia1(j)) - end do - enddo - else - k=0 - - do i=irw,lrw - do j=a%ia2(i),a%ia2(i+1)-1 - k = k + 1 - ia(nzin_+k) = i - ja(nzin_+k) = a%ia1(j) - val(nzin_+k) = a%aspk(j) - end do - enddo - end if -!!$ if (nz /= k) then -!!$ write(0,*) 'csr getrow Size mismatch ',nz,k +!!$ subroutine csr_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) +!!$ +!!$ use psb_sort_mod +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(in) :: a +!!$ integer :: irw +!!$ integer, intent(out) :: nz +!!$ integer, allocatable, intent(inout) :: ia(:), ja(:) +!!$ real(psb_dpk_), allocatable, intent(inout) :: val(:) +!!$ integer :: nzin +!!$ logical, intent(in) :: append +!!$ integer :: lrw,info +!!$ integer, optional :: iren(:) +!!$ +!!$ integer :: idx,i,j, k, nr, row_idx, nzin_ +!!$ integer, allocatable :: indices(:) +!!$ +!!$ if (append) then +!!$ nzin_ = nzin +!!$ else +!!$ nzin_ = 0 +!!$ endif +!!$ +!!$ if (a%pl(1) /= 0) then +!!$ +!!$ nr = lrw - irw + 1 +!!$ allocate(indices(nr)) +!!$ nz = 0 +!!$ do i=1,nr +!!$ indices(i)=a%pl(irw+i-1) +!!$ nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i)) +!!$ end do +!!$ +!!$ call psb_ensure_size(nzin_+nz,ia,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,ja,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,val,info) +!!$ if (info /= 0) return +!!$ +!!$ k=0 +!!$ if(present(iren)) then +!!$ do i=1,nr +!!$ row_idx=indices(i) +!!$ do j=a%ia2(row_idx),a%ia2(row_idx+1)-1 +!!$ k = k + 1 +!!$ val(nzin_+k) = a%aspk(j) +!!$ ia(nzin_+k) = iren(row_idx) +!!$ ja(nzin_+k) = iren(a%ia1(j)) +!!$ end do +!!$ end do +!!$ else +!!$ do i=1,nr +!!$ row_idx=indices(i) +!!$ do j=a%ia2(row_idx),a%ia2(row_idx+1)-1 +!!$ k = k + 1 +!!$ val(nzin_+k) = a%aspk(j) +!!$ ia(nzin_+k) = row_idx +!!$ ja(nzin_+k) = a%ia1(j) +!!$ end do +!!$ end do +!!$ end if +!!$ +!!$ else +!!$ +!!$ idx = irw +!!$ +!!$ if (idx<0) then +!!$ write(0,*) ' spgtrow Error : idx no good ',idx +!!$ info = 2 +!!$ return +!!$ end if +!!$ nr = lrw - irw + 1 +!!$ nz = a%ia2(idx+nr) - a%ia2(idx) +!!$ +!!$ call psb_ensure_size(nzin_+nz,ia,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,ja,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,val,info) +!!$ if (info /= 0) return +!!$ +!!$ +!!$ if (present(iren)) then +!!$ k=0 +!!$ do i=irw,lrw +!!$ do j=a%ia2(i),a%ia2(i+1)-1 +!!$ k = k + 1 +!!$ val(nzin_+k) = a%aspk(j) +!!$ ia(nzin_+k) = iren(i) +!!$ ja(nzin_+k) = iren(a%ia1(j)) +!!$ end do +!!$ enddo +!!$ else +!!$ k=0 +!!$ +!!$ do i=irw,lrw +!!$ do j=a%ia2(i),a%ia2(i+1)-1 +!!$ k = k + 1 +!!$ ia(nzin_+k) = i +!!$ ja(nzin_+k) = a%ia1(j) +!!$ val(nzin_+k) = a%aspk(j) +!!$ end do +!!$ enddo +!!$ end if +! !$ if (nz /= k) then +! !$ write(0,*) 'csr getrow Size mismatch ',nz,k +! !$ endif +!!$ if (a%pr(1) /= 0) then +!!$ write(0,*) 'Feeling lazy today, Right Permutation will have to wait' !!$ endif - if (a%pr(1) /= 0) then - write(0,*) 'Feeling lazy today, Right Permutation will have to wait' - endif - - endif - - end subroutine csr_dspgtrow - - subroutine coo_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) - - use psb_sort_mod - use psb_spmat_type - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - integer :: irw - integer, intent(out) :: nz - integer, allocatable, intent(inout) :: ia(:), ja(:) - real(psb_dpk_), allocatable, intent(inout) :: val(:) - integer :: nzin - logical, intent(in) :: append - integer :: lrw,info - integer, optional :: iren(:) - integer :: nzin_, nza, idx,ip,jp,i,k, nzt - integer :: debug_level, debug_unit - character(len=20) :: name='coo_dspgtrow' - - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - nza = a%infoa(psb_nnz_) - if (a%pl(1) /= 0) then - write(debug_unit,*) 'Fatal error in SPGTROW: do not feed a permuted mat so far!' - idx = -1 - else - idx = irw - endif - if (idx<0) then - write(debug_unit,*) ' spgtrow Error : idx no good ',idx - info = 2 - return - end if - - if (append) then - nzin_ = nzin - else - nzin_ = 0 - endif - - if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then - ! In this case we can do a binary search. - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name), ': srtdcoo ' - do - ip = psb_ibsrch(irw,nza,a%ia1) - if (ip /= -1) exit - irw = irw + 1 - if (irw > lrw) then - write(debug_unit,*) trim(name),& - & 'Warning : did not find any rows. Is this an error? ',& - & irw,lrw,idx - exit - end if - end do - - if (ip /= -1) then - ! expand [ip,jp] to contain all row entries. - do - if (ip < 2) exit - if (a%ia1(ip-1) == irw) then - ip = ip -1 - else - exit - end if - end do - - end if - - do - jp = psb_ibsrch(lrw,nza,a%ia1) - 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%ia1(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. - nz = jp - ip +1 - - call psb_ensure_size(nzin_+nz,ia,info) - if (info==0) call psb_ensure_size(nzin_+nz,ja,info) - if (info==0) call psb_ensure_size(nzin_+nz,val,info) - if (info /= 0) return - - if (present(iren)) then - do i=ip,jp - nzin_ = nzin_ + 1 - val(nzin_) = a%aspk(i) - ia(nzin_) = iren(a%ia1(i)) - ja(nzin_) = iren(a%ia2(i)) - enddo - else - do i=ip,jp - nzin_ = nzin_ + 1 - val(nzin_) = a%aspk(i) - ia(nzin_) = a%ia1(i) - ja(nzin_) = a%ia2(i) - 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%m,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%infoa(psb_nnz_) - if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) 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%aspk(i) - ia(nzin_+k) = iren(a%ia1(i)) - ja(nzin_+k) = iren(a%ia2(i)) - endif - enddo - else - k = 0 - do i=1,a%infoa(psb_nnz_) - if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) 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%aspk(i) - ia(nzin_+k) = (a%ia1(i)) - ja(nzin_+k) = (a%ia2(i)) - endif - enddo - nzin_=nzin_+k - end if - nz = k - end if - - end subroutine coo_dspgtrow - - - subroutine jad_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) - - use psb_sort_mod - use psb_spmat_type - use psb_const_mod - - implicit none - - type(psb_dspmat_type), intent(in), target :: a - integer :: irw - integer, intent(out) :: nz - integer, allocatable, intent(inout) :: ia(:), ja(:) - real(psb_dpk_), allocatable, intent(inout) :: val(:) - integer :: nzin - logical, intent(in) :: append - integer, optional :: iren(:) - integer :: lrw,info - - integer, pointer :: ia1(:), ia2(:), ia3(:),& - & ja_(:), ka_(:), indices(:), blks(:) - integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng, nzin_,& - & i,j,k,nr - - - png = a%ia2(1) ! points to the number of blocks - pia = a%ia2(2) ! points to the beginning of ia(3,png) - pja = a%ia2(3) ! points to the beginning of ja(:) - - ng = a%ia2(png) ! the number of blocks - ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk - ka_ => a%ia1(:) ! the array containing the column indices - ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block - ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column - ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column - - if (append) then - nzin_ = nzin - else - nzin_ = 0 - endif - - if (a%pl(1) /= 0) then - - nr = lrw - irw + 1 - allocate(indices(nr),blks(nr)) - nz = 0 - - do i=1,nr - indices(i)=a%pl(irw+i-1) - j=0 - blkfnd: do - j=j+1 - if(ia1(j) == indices(i)) then - blks(i)=j - nz=nz+ia3(j)-ia2(j) - ipx = ia1(j) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j)+rb - nz = nz+ja_(row+1)-ja_(row) - exit blkfnd - else if(ia1(j) > indices(i)) then - blks(i)=j-1 - nz=nz+ia3(j-1)-ia2(j-1) - ipx = ia1(j-1) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j-1)+rb - nz = nz+ja_(row+1)-ja_(row) - exit blkfnd - end if - end do blkfnd - end do - - - call psb_ensure_size(nzin_+nz,ia,info) - if (info==0) call psb_ensure_size(nzin_+nz,ja,info) - if (info==0) call psb_ensure_size(nzin_+nz,val,info) - if (info /= 0) return - - k=0 - ! cycle over rows - do i=1,nr - - ! find which block the row belongs to - blk = blks(i) - - ! extract first part of the row from the jad block - ipx = ia1(blk) ! the first row index of the block - k_pt= ia2(blk) ! the pointer to the beginning of a column in ja - rb = indices(i)-ipx ! the row offset within the block - npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block - - if(present(iren))then - do col = ia2(blk), ia3(blk)-1 - k=k+1 - val(nzin_+k) = a%aspk(ja_(col)+rb) - ia(nzin_+k) = iren(irw+i-1) - ja(nzin_+k) = iren(ka_(ja_(col)+rb)) - end do - else - do col = ia2(blk), ia3(blk)-1 - k=k+1 - val(nzin_+k) = a%aspk(ja_(col)+rb) - ia(nzin_+k) = irw+i-1 - ja(nzin_+k) = ka_(ja_(col)+rb) - end do - end if - ! extract second part of the row from the csr tail - row=ia3(blk)+rb - if(present(iren))then - do j=ja_(row), ja_(row+1)-1 - k=k+1 - val(nzin_+k) = a%aspk(j) - ia(nzin_+k) = iren(irw+i-1) - ja(nzin_+k) = iren(ka_(j)) - end do - else - do j=ja_(row), ja_(row+1)-1 - k=k+1 - val(nzin_+k) = a%aspk(j) - ia(nzin_+k) = irw+i-1 - ja(nzin_+k) = ka_(j) - end do - end if - end do - - else - ! There might be some problems - info=134 - end if - - end subroutine jad_dspgtrow +!!$ +!!$ endif +!!$ +!!$ end subroutine csr_dspgtrow + +!!$ subroutine coo_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) +!!$ +!!$ use psb_sort_mod +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(in) :: a +!!$ integer :: irw +!!$ integer, intent(out) :: nz +!!$ integer, allocatable, intent(inout) :: ia(:), ja(:) +!!$ real(psb_dpk_), allocatable, intent(inout) :: val(:) +!!$ integer :: nzin +!!$ logical, intent(in) :: append +!!$ integer :: lrw,info +!!$ integer, optional :: iren(:) +!!$ integer :: nzin_, nza, idx,ip,jp,i,k, nzt +!!$ integer :: debug_level, debug_unit +!!$ character(len=20) :: name='coo_dspgtrow' +!!$ +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ nza = a%infoa(psb_nnz_) +!!$ if (a%pl(1) /= 0) then +!!$ write(debug_unit,*) 'Fatal error in SPGTROW: do not feed a permuted mat so far!' +!!$ idx = -1 +!!$ else +!!$ idx = irw +!!$ endif +!!$ if (idx<0) then +!!$ write(debug_unit,*) ' spgtrow Error : idx no good ',idx +!!$ info = 2 +!!$ return +!!$ end if +!!$ +!!$ if (append) then +!!$ nzin_ = nzin +!!$ else +!!$ nzin_ = 0 +!!$ endif +!!$ +!!$ if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then +!!$ ! In this case we can do a binary search. +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name), ': srtdcoo ' +!!$ do +!!$ ip = psb_ibsrch(irw,nza,a%ia1) +!!$ if (ip /= -1) exit +!!$ irw = irw + 1 +!!$ if (irw > lrw) then +!!$ write(debug_unit,*) trim(name),& +!!$ & 'Warning : did not find any rows. Is this an error? ',& +!!$ & irw,lrw,idx +!!$ exit +!!$ end if +!!$ end do +!!$ +!!$ if (ip /= -1) then +!!$ ! expand [ip,jp] to contain all row entries. +!!$ do +!!$ if (ip < 2) exit +!!$ if (a%ia1(ip-1) == irw) then +!!$ ip = ip -1 +!!$ else +!!$ exit +!!$ end if +!!$ end do +!!$ +!!$ end if +!!$ +!!$ do +!!$ jp = psb_ibsrch(lrw,nza,a%ia1) +!!$ 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%ia1(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. +!!$ nz = jp - ip +1 +!!$ +!!$ call psb_ensure_size(nzin_+nz,ia,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,ja,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,val,info) +!!$ if (info /= 0) return +!!$ +!!$ if (present(iren)) then +!!$ do i=ip,jp +!!$ nzin_ = nzin_ + 1 +!!$ val(nzin_) = a%aspk(i) +!!$ ia(nzin_) = iren(a%ia1(i)) +!!$ ja(nzin_) = iren(a%ia2(i)) +!!$ enddo +!!$ else +!!$ do i=ip,jp +!!$ nzin_ = nzin_ + 1 +!!$ val(nzin_) = a%aspk(i) +!!$ ia(nzin_) = a%ia1(i) +!!$ ja(nzin_) = a%ia2(i) +!!$ 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%m,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%infoa(psb_nnz_) +!!$ if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) 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%aspk(i) +!!$ ia(nzin_+k) = iren(a%ia1(i)) +!!$ ja(nzin_+k) = iren(a%ia2(i)) +!!$ endif +!!$ enddo +!!$ else +!!$ k = 0 +!!$ do i=1,a%infoa(psb_nnz_) +!!$ if ((a%ia1(i)>=irw).and.(a%ia1(i)<=lrw)) 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%aspk(i) +!!$ ia(nzin_+k) = (a%ia1(i)) +!!$ ja(nzin_+k) = (a%ia2(i)) +!!$ endif +!!$ enddo +!!$ nzin_=nzin_+k +!!$ end if +!!$ nz = k +!!$ end if +!!$ +!!$ end subroutine coo_dspgtrow +!!$ +!!$ +!!$ subroutine jad_dspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) +!!$ +!!$ use psb_sort_mod +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(in), target :: a +!!$ integer :: irw +!!$ integer, intent(out) :: nz +!!$ integer, allocatable, intent(inout) :: ia(:), ja(:) +!!$ real(psb_dpk_), allocatable, intent(inout) :: val(:) +!!$ integer :: nzin +!!$ logical, intent(in) :: append +!!$ integer, optional :: iren(:) +!!$ integer :: lrw,info +!!$ +!!$ integer, pointer :: ia1(:), ia2(:), ia3(:),& +!!$ & ja_(:), ka_(:), indices(:), blks(:) +!!$ integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng, nzin_,& +!!$ & i,j,k,nr +!!$ +!!$ +!!$ png = a%ia2(1) ! points to the number of blocks +!!$ pia = a%ia2(2) ! points to the beginning of ia(3,png) +!!$ pja = a%ia2(3) ! points to the beginning of ja(:) +!!$ +!!$ ng = a%ia2(png) ! the number of blocks +!!$ ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk +!!$ ka_ => a%ia1(:) ! the array containing the column indices +!!$ ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index of each block +!!$ ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. in ja to the first jad column +!!$ ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. in ja to the first csr column +!!$ +!!$ if (append) then +!!$ nzin_ = nzin +!!$ else +!!$ nzin_ = 0 +!!$ endif +!!$ +!!$ if (a%pl(1) /= 0) then +!!$ +!!$ nr = lrw - irw + 1 +!!$ allocate(indices(nr),blks(nr)) +!!$ nz = 0 +!!$ +!!$ do i=1,nr +!!$ indices(i)=a%pl(irw+i-1) +!!$ j=0 +!!$ blkfnd: do +!!$ j=j+1 +!!$ if(ia1(j) == indices(i)) then +!!$ blks(i)=j +!!$ nz=nz+ia3(j)-ia2(j) +!!$ ipx = ia1(j) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j)+rb +!!$ nz = nz+ja_(row+1)-ja_(row) +!!$ exit blkfnd +!!$ else if(ia1(j) > indices(i)) then +!!$ blks(i)=j-1 +!!$ nz=nz+ia3(j-1)-ia2(j-1) +!!$ ipx = ia1(j-1) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j-1)+rb +!!$ nz = nz+ja_(row+1)-ja_(row) +!!$ exit blkfnd +!!$ end if +!!$ end do blkfnd +!!$ end do +!!$ +!!$ +!!$ call psb_ensure_size(nzin_+nz,ia,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,ja,info) +!!$ if (info==0) call psb_ensure_size(nzin_+nz,val,info) +!!$ if (info /= 0) return +!!$ +!!$ k=0 +!!$ ! cycle over rows +!!$ do i=1,nr +!!$ +!!$ ! find which block the row belongs to +!!$ blk = blks(i) +!!$ +!!$ ! extract first part of the row from the jad block +!!$ ipx = ia1(blk) ! the first row index of the block +!!$ k_pt= ia2(blk) ! the pointer to the beginning of a column in ja +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block +!!$ +!!$ if(present(iren))then +!!$ do col = ia2(blk), ia3(blk)-1 +!!$ k=k+1 +!!$ val(nzin_+k) = a%aspk(ja_(col)+rb) +!!$ ia(nzin_+k) = iren(irw+i-1) +!!$ ja(nzin_+k) = iren(ka_(ja_(col)+rb)) +!!$ end do +!!$ else +!!$ do col = ia2(blk), ia3(blk)-1 +!!$ k=k+1 +!!$ val(nzin_+k) = a%aspk(ja_(col)+rb) +!!$ ia(nzin_+k) = irw+i-1 +!!$ ja(nzin_+k) = ka_(ja_(col)+rb) +!!$ end do +!!$ end if +!!$ ! extract second part of the row from the csr tail +!!$ row=ia3(blk)+rb +!!$ if(present(iren))then +!!$ do j=ja_(row), ja_(row+1)-1 +!!$ k=k+1 +!!$ val(nzin_+k) = a%aspk(j) +!!$ ia(nzin_+k) = iren(irw+i-1) +!!$ ja(nzin_+k) = iren(ka_(j)) +!!$ end do +!!$ else +!!$ do j=ja_(row), ja_(row+1)-1 +!!$ k=k+1 +!!$ val(nzin_+k) = a%aspk(j) +!!$ ia(nzin_+k) = irw+i-1 +!!$ ja(nzin_+k) = ka_(j) +!!$ end do +!!$ end if +!!$ end do +!!$ +!!$ else +!!$ ! There might be some problems +!!$ info=134 +!!$ end if +!!$ +!!$ end subroutine jad_dspgtrow subroutine csr_cspgtrow(irw,a,nz,ia,ja,val,nzin,append,lrw,info,iren) diff --git a/base/serial/psb_regen_mod.f90 b/base/serial/psb_regen_mod.f90 index 76f0a06a..ffb9e5eb 100644 --- a/base/serial/psb_regen_mod.f90 +++ b/base/serial/psb_regen_mod.f90 @@ -32,13 +32,13 @@ module psb_regen_mod interface csr_regen - module procedure csr_ssp_regen, csr_dsp_regen, csr_csp_regen, csr_zsp_regen + module procedure csr_ssp_regen, csr_csp_regen, csr_zsp_regen end interface interface coo_regen - module procedure coo_ssp_regen, coo_dsp_regen, coo_csp_regen, coo_zsp_regen + module procedure coo_ssp_regen, coo_csp_regen, coo_zsp_regen end interface interface jad_regen - module procedure jad_ssp_regen, jad_dsp_regen, jad_csp_regen, jad_zsp_regen + module procedure jad_ssp_regen, jad_csp_regen, jad_zsp_regen end interface contains @@ -360,323 +360,323 @@ contains end subroutine jad_ssp_regen - subroutine csr_dsp_regen(a,info) - - use psb_spmat_type - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer :: info - - integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt - real(psb_dpk_), allocatable :: work(:) - integer :: err_act - character(len=20) :: name - integer :: debug_level, debug_unit - - - name='psb_spcnv' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - - ! - ! dupl_ and upd_ fields should not be changed. - ! - select case(psb_sp_getifld(psb_upd_,a,info)) - - case(psb_upd_perm_) - - allocate(work(size(a%aspk)+1000),stat=info) - if (info /= 0) then - info=2040 - call psb_errpush(info,name) - goto 9999 - end if - - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' - ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) - ip2 = a%ia2(ip1+psb_ip2_) - nnz = a%ia2(ip1+psb_nnz_) - iflag = a%ia2(ip1+psb_iflag_) - ichk = a%ia2(ip1+psb_ichk_) - nnzt = a%ia2(ip1+psb_nnzt_) - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),'Regeneration start: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info - - if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then - info = 8889 - write(debug_unit,*) trim(name),'Regeneration start error: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk - call psb_errpush(info,name) - goto 9999 - endif - do i= 1, nnz - work(i) = dzero - enddo - select case(iflag) - case(psb_dupl_ovwrt_,psb_dupl_err_) - do i=1, nnz - work(a%ia2(ip2+i-1)) = a%aspk(i) - enddo - case(psb_dupl_add_) - do i=1, nnz - work(a%ia2(ip2+i-1)) = a%aspk(i) + work(a%ia2(ip2+i-1)) - enddo - case default - info = 8887 - call psb_errpush(info,name) - goto 9999 - end select - - do i=1, nnz - a%aspk(i) = work(i) - enddo - - - case(psb_upd_srch_) - ! Nothing to be done here. - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),& - & 'Going through on regeneration with psb_upd_srch_' - case default - ! Wrong value - info = 8888 - call psb_errpush(info,name) - goto 9999 - - end select - - 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 csr_dsp_regen - - subroutine coo_dsp_regen(a,info) - - use psb_spmat_type - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer :: info - - integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt - real(psb_dpk_), allocatable :: work(:) - integer :: err_act - character(len=20) :: name - integer :: debug_level, debug_unit - - - name='psb_spcnv' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - - ! - ! dupl_ and upd_ fields should not be changed. - ! - select case(psb_sp_getifld(psb_upd_,a,info)) - - case(psb_upd_perm_) - - allocate(work(size(a%aspk)+1000),stat=info) - if (info /= 0) then - info=2040 - call psb_errpush(info,name) - goto 9999 - end if - - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' - ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) - ip2 = a%ia2(ip1+psb_ip2_) - nnz = a%ia2(ip1+psb_nnz_) - iflag = a%ia2(ip1+psb_iflag_) - ichk = a%ia2(ip1+psb_ichk_) - nnzt = a%ia2(ip1+psb_nnzt_) - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),'Regeneration start: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info - - if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then - info = 8889 - write(debug_unit,*) trim(name),'Regeneration start error: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk - call psb_errpush(info,name) - goto 9999 - endif - do i= 1, nnz - work(i) = dzero - enddo - select case(iflag) - case(psb_dupl_ovwrt_,psb_dupl_err_) - do i=1, nnz - work(a%ia2(ip2+i-1)) = a%aspk(i) - enddo - case(psb_dupl_add_) - do i=1, nnz - work(a%ia2(ip2+i-1)) = a%aspk(i) + work(a%ia2(ip2+i-1)) - enddo - case default - info = 8887 - call psb_errpush(info,name) - goto 9999 - end select - - do i=1, nnz - a%aspk(i) = work(i) - enddo - - - case(psb_upd_srch_) - ! Nothing to be done here. - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & 'Going through on regeneration with psb_upd_srch_' - case default - ! Wrong value - info = 8888 - call psb_errpush(info,name) - goto 9999 - - end select - - 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 coo_dsp_regen - - subroutine jad_dsp_regen(a,info) - - use psb_spmat_type - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer :: info - - integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt - real(psb_dpk_), allocatable :: work(:) - integer :: err_act - character(len=20) :: name - integer :: debug_level, debug_unit - - name='psb_spcnv' - info = 0 - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - - ! - ! dupl_ and upd_ fields should not be changed. - ! - select case(psb_sp_getifld(psb_upd_,a,info)) - - case(psb_upd_perm_) - - allocate(work(size(a%aspk)+1000),stat=info) - if (info /= 0) then - info=2040 - call psb_errpush(info,name) - goto 9999 - end if - - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' - ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) - ip2 = a%ia1(ip1+psb_ip2_) - nnz = a%ia1(ip1+psb_nnz_) - iflag = a%ia1(ip1+psb_iflag_) - ichk = a%ia1(ip1+psb_ichk_) - nnzt = a%ia1(ip1+psb_nnzt_) - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),'Regeneration start: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info - - if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then - info = 8889 - write(debug_unit,*) trim(name),'Regeneration start error: ',& - & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk - call psb_errpush(info,name) - goto 9999 - endif - do i= 1, nnz - work(i) = dzero - enddo - select case(iflag) - case(psb_dupl_ovwrt_,psb_dupl_err_) - do i=1, nnz - work(a%ia1(ip2+i-1)) = a%aspk(i) - enddo - case(psb_dupl_add_) - do i=1, nnz - work(a%ia1(ip2+i-1)) = a%aspk(i) + work(a%ia1(ip2+i-1)) - enddo - case default - info = 8887 - call psb_errpush(info,name) - goto 9999 - end select - - do i=1, nnz - a%aspk(i) = work(i) - enddo - - - case(psb_upd_srch_) - ! Nothing to be done here. - if (debug_level >= psb_debug_serial_)& - & write(debug_unit,*) trim(name),& - & 'Going through on regeneration with psb_upd_srch_' - case default - ! Wrong value - info = 8888 - call psb_errpush(info,name) - goto 9999 - - end select - - 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 jad_dsp_regen - +!!$ subroutine csr_dsp_regen(a,info) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer :: info +!!$ +!!$ integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt +!!$ real(psb_dpk_), allocatable :: work(:) +!!$ integer :: err_act +!!$ character(len=20) :: name +!!$ integer :: debug_level, debug_unit +!!$ +!!$ +!!$ name='psb_spcnv' +!!$ info = 0 +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ +!!$ ! +!!$ ! dupl_ and upd_ fields should not be changed. +!!$ ! +!!$ select case(psb_sp_getifld(psb_upd_,a,info)) +!!$ +!!$ case(psb_upd_perm_) +!!$ +!!$ allocate(work(size(a%aspk)+1000),stat=info) +!!$ if (info /= 0) then +!!$ info=2040 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end if +!!$ +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' +!!$ ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) +!!$ ip2 = a%ia2(ip1+psb_ip2_) +!!$ nnz = a%ia2(ip1+psb_nnz_) +!!$ iflag = a%ia2(ip1+psb_iflag_) +!!$ ichk = a%ia2(ip1+psb_ichk_) +!!$ nnzt = a%ia2(ip1+psb_nnzt_) +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),'Regeneration start: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info +!!$ +!!$ if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then +!!$ info = 8889 +!!$ write(debug_unit,*) trim(name),'Regeneration start error: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ endif +!!$ do i= 1, nnz +!!$ work(i) = dzero +!!$ enddo +!!$ select case(iflag) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ do i=1, nnz +!!$ work(a%ia2(ip2+i-1)) = a%aspk(i) +!!$ enddo +!!$ case(psb_dupl_add_) +!!$ do i=1, nnz +!!$ work(a%ia2(ip2+i-1)) = a%aspk(i) + work(a%ia2(ip2+i-1)) +!!$ enddo +!!$ case default +!!$ info = 8887 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, nnz +!!$ a%aspk(i) = work(i) +!!$ enddo +!!$ +!!$ +!!$ case(psb_upd_srch_) +!!$ ! Nothing to be done here. +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),& +!!$ & 'Going through on regeneration with psb_upd_srch_' +!!$ case default +!!$ ! Wrong value +!!$ info = 8888 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ +!!$ end select +!!$ +!!$ 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 csr_dsp_regen +!!$ +!!$ subroutine coo_dsp_regen(a,info) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer :: info +!!$ +!!$ integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt +!!$ real(psb_dpk_), allocatable :: work(:) +!!$ integer :: err_act +!!$ character(len=20) :: name +!!$ integer :: debug_level, debug_unit +!!$ +!!$ +!!$ name='psb_spcnv' +!!$ info = 0 +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ +!!$ ! +!!$ ! dupl_ and upd_ fields should not be changed. +!!$ ! +!!$ select case(psb_sp_getifld(psb_upd_,a,info)) +!!$ +!!$ case(psb_upd_perm_) +!!$ +!!$ allocate(work(size(a%aspk)+1000),stat=info) +!!$ if (info /= 0) then +!!$ info=2040 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end if +!!$ +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' +!!$ ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) +!!$ ip2 = a%ia2(ip1+psb_ip2_) +!!$ nnz = a%ia2(ip1+psb_nnz_) +!!$ iflag = a%ia2(ip1+psb_iflag_) +!!$ ichk = a%ia2(ip1+psb_ichk_) +!!$ nnzt = a%ia2(ip1+psb_nnzt_) +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),'Regeneration start: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info +!!$ +!!$ if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then +!!$ info = 8889 +!!$ write(debug_unit,*) trim(name),'Regeneration start error: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ endif +!!$ do i= 1, nnz +!!$ work(i) = dzero +!!$ enddo +!!$ select case(iflag) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ do i=1, nnz +!!$ work(a%ia2(ip2+i-1)) = a%aspk(i) +!!$ enddo +!!$ case(psb_dupl_add_) +!!$ do i=1, nnz +!!$ work(a%ia2(ip2+i-1)) = a%aspk(i) + work(a%ia2(ip2+i-1)) +!!$ enddo +!!$ case default +!!$ info = 8887 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, nnz +!!$ a%aspk(i) = work(i) +!!$ enddo +!!$ +!!$ +!!$ case(psb_upd_srch_) +!!$ ! Nothing to be done here. +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & 'Going through on regeneration with psb_upd_srch_' +!!$ case default +!!$ ! Wrong value +!!$ info = 8888 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ +!!$ end select +!!$ +!!$ 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 coo_dsp_regen +!!$ +!!$ subroutine jad_dsp_regen(a,info) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer :: info +!!$ +!!$ integer :: i, ip1,ip2,nnz,iflag,ichk,nnzt +!!$ real(psb_dpk_), allocatable :: work(:) +!!$ integer :: err_act +!!$ character(len=20) :: name +!!$ integer :: debug_level, debug_unit +!!$ +!!$ name='psb_spcnv' +!!$ info = 0 +!!$ call psb_erractionsave(err_act) +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ +!!$ ! +!!$ ! dupl_ and upd_ fields should not be changed. +!!$ ! +!!$ select case(psb_sp_getifld(psb_upd_,a,info)) +!!$ +!!$ case(psb_upd_perm_) +!!$ +!!$ allocate(work(size(a%aspk)+1000),stat=info) +!!$ if (info /= 0) then +!!$ info=2040 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end if +!!$ +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),'Regeneration with psb_upd_perm_' +!!$ ip1 = psb_sp_getifld(psb_upd_pnt_,a,info) +!!$ ip2 = a%ia1(ip1+psb_ip2_) +!!$ nnz = a%ia1(ip1+psb_nnz_) +!!$ iflag = a%ia1(ip1+psb_iflag_) +!!$ ichk = a%ia1(ip1+psb_ichk_) +!!$ nnzt = a%ia1(ip1+psb_nnzt_) +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),'Regeneration start: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,info +!!$ +!!$ if ((ichk/=nnzt+iflag).or.(nnz/=nnzt)) then +!!$ info = 8889 +!!$ write(debug_unit,*) trim(name),'Regeneration start error: ',& +!!$ & a%infoa(psb_upd_),psb_upd_perm_,nnz,nnzt ,iflag,ichk +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ endif +!!$ do i= 1, nnz +!!$ work(i) = dzero +!!$ enddo +!!$ select case(iflag) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ do i=1, nnz +!!$ work(a%ia1(ip2+i-1)) = a%aspk(i) +!!$ enddo +!!$ case(psb_dupl_add_) +!!$ do i=1, nnz +!!$ work(a%ia1(ip2+i-1)) = a%aspk(i) + work(a%ia1(ip2+i-1)) +!!$ enddo +!!$ case default +!!$ info = 8887 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, nnz +!!$ a%aspk(i) = work(i) +!!$ enddo +!!$ +!!$ +!!$ case(psb_upd_srch_) +!!$ ! Nothing to be done here. +!!$ if (debug_level >= psb_debug_serial_)& +!!$ & write(debug_unit,*) trim(name),& +!!$ & 'Going through on regeneration with psb_upd_srch_' +!!$ case default +!!$ ! Wrong value +!!$ info = 8888 +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ +!!$ end select +!!$ +!!$ 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 jad_dsp_regen +!!$ subroutine csr_csp_regen(a,info) diff --git a/base/serial/psb_update_mod.f90 b/base/serial/psb_update_mod.f90 index 709ce25a..c028f5a3 100644 --- a/base/serial/psb_update_mod.f90 +++ b/base/serial/psb_update_mod.f90 @@ -32,19 +32,19 @@ module psb_update_mod interface psb_srch_upd - module procedure psb_s_srch_upd, psb_d_srch_upd, psb_c_srch_upd, psb_z_srch_upd + module procedure psb_s_srch_upd, psb_c_srch_upd, psb_z_srch_upd end interface interface coo_srch_upd - module procedure s_coo_srch_upd, d_coo_srch_upd, c_coo_srch_upd, z_coo_srch_upd + module procedure s_coo_srch_upd, c_coo_srch_upd, z_coo_srch_upd end interface interface csr_srch_upd - module procedure s_csr_srch_upd, d_csr_srch_upd, c_csr_srch_upd, z_csr_srch_upd + module procedure s_csr_srch_upd, c_csr_srch_upd, z_csr_srch_upd end interface interface jad_srch_upd - module procedure s_jad_srch_upd, d_jad_srch_upd, c_jad_srch_upd, z_jad_srch_upd + module procedure s_jad_srch_upd, c_jad_srch_upd, z_jad_srch_upd end interface contains @@ -95,50 +95,50 @@ contains end subroutine psb_s_srch_upd - subroutine psb_d_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - - use psb_spmat_type - use psb_const_mod - use psb_realloc_mod - use psb_string_mod - use psb_serial_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl - integer, intent(in) :: ia(:),ja(:) - integer, intent(inout) :: nza - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: ng,gtl(:) - - info = 0 - - if (present(gtl)) then - if (.not.present(ng)) then - info = -1 - return - endif - end if - select case(psb_tolower(a%fida)) - case ('csr') - call csr_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - case ('coo') - call coo_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - case ('jad') - call jad_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - - case default - - info = -9 - - end select - - end subroutine psb_d_srch_upd +!!$ subroutine psb_d_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_realloc_mod +!!$ use psb_string_mod +!!$ use psb_serial_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ integer, intent(inout) :: nza +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: ng,gtl(:) +!!$ +!!$ info = 0 +!!$ +!!$ if (present(gtl)) then +!!$ if (.not.present(ng)) then +!!$ info = -1 +!!$ return +!!$ endif +!!$ end if +!!$ select case(psb_tolower(a%fida)) +!!$ case ('csr') +!!$ call csr_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ case ('coo') +!!$ call coo_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ case ('jad') +!!$ call jad_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ +!!$ case default +!!$ +!!$ info = -9 +!!$ +!!$ end select +!!$ +!!$ end subroutine psb_d_srch_upd subroutine psb_c_srch_upd(nz,ia,ja,val,nza,a,& & imin,imax,jmin,jmax,nzl,info,gtl,ng) @@ -915,658 +915,658 @@ contains end subroutine s_jad_srch_upd - subroutine d_csr_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - - use psb_spmat_type - use psb_const_mod - use psb_realloc_mod - use psb_string_mod - use psb_serial_mod - use psb_error_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl - integer, intent(in) :: ia(:),ja(:) - integer, intent(inout) :: nza - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: ng,gtl(:) - - integer :: debug_level, debug_unit - character(len=20) :: name='d_csr_srch_upd' - integer :: i,ir,ic, ilr, ilc, ip, & - & i1,i2,nc,lb,ub,m,dupl - - info = 0 - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - dupl = psb_sp_getifld(psb_dupl_,a,info) - - if (present(gtl)) then - if (.not.present(ng)) then - info = -1 - return - endif - - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then - i1 = a%ia2(ir) - i2 = a%ia2(ir+1) - nc=i2-i1 - - ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) - if (ip>0) then - a%aspk(i1+ip-1) = val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',a%ia1(i1:i2-1) - info = i - return - end if - - else - - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - end if - end do - - case(psb_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then - i1 = a%ia2(ir) - i2 = a%ia2(ir+1) - nc = i2-i1 - ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) - if (ip>0) then - a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',a%ia1(i1:i2-1) - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - - else - - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - - if ((ir > 0).and.(ir <= a%m)) then - - i1 = a%ia2(ir) - i2 = a%ia2(ir+1) - nc=i2-i1 - - ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) - if (ip>0) then - a%aspk(i1+ip-1) = val(i) - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Was searching ',ic,' in: ',i1,i2,& - & ' : ',a%ia1(i1:i2-1) - info = i - return - end if - - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - - end do - - case(psb_dupl_add_) - ! Add - ilr = -1 - ilc = -1 - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then - i1 = a%ia2(ir) - i2 = a%ia2(ir+1) - nc = i2-i1 - ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) - if (ip>0) then - a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) - else - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - end if - - end subroutine d_csr_srch_upd - - subroutine d_coo_srch_upd(nz,ia,ja,val,nza,a,& - & imin,imax,jmin,jmax,nzl,info,gtl,ng) - - use psb_spmat_type - use psb_const_mod - use psb_realloc_mod - use psb_string_mod - use psb_serial_mod - implicit none - - type(psb_dspmat_type), intent(inout) :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl - integer, intent(in) :: ia(:),ja(:) - integer, intent(inout) :: nza - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: ng,gtl(:) - integer :: i,ir,ic, ilr, ilc, ip, & - & i1,i2,nc,nnz,dupl - integer :: debug_level, debug_unit - character(len=20) :: name='d_coo_srch_upd' - - info = 0 - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - dupl = psb_sp_getifld(psb_dupl_,a,info) - - if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then - info = -4 - return - end if - - ilr = -1 - ilc = -1 - nnz = psb_sp_getifld(psb_nnz_,a,info) - - - if (present(gtl)) then - if (.not.present(ng)) then - info = -1 - return - endif - - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - if ((ir > 0).and.(ir <= a%m)) then - ic = gtl(ic) - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nnz,a%ia1) - i2 = i1 - do - if (i2+1 > nnz) exit - if (a%ia1(i2+1) /= a%ia1(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (a%ia1(i1-1) /= a%ia1(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,a%ia2(i1:i2)) - if (ip>0) then - a%aspk(i1+ip-1) = val(i) - else - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - endif - end if - end do - case(psb_dupl_add_) - ! Add - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - ic = gtl(ic) - if ((ir > 0).and.(ir <= a%m)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nnz,a%ia1) - i2 = i1 - do - if (i2+1 > nnz) exit - if (a%ia1(i2+1) /= a%ia1(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (a%ia1(i1-1) /= a%ia1(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,a%ia2(i1:i2)) - if (ip>0) then - a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) - else - info = i - return - end if - else - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Discarding row that does not belong to us.' - end if - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - - else - - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nnz,a%ia1) - i2 = i1 - do - if (i2+1 > nnz) exit - if (a%ia1(i2+1) /= a%ia1(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (a%ia1(i1-1) /= a%ia1(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,a%ia2(i1:i2)) - if (ip>0) then - a%aspk(i1+ip-1) = val(i) - else - info = i - return - end if - end if - end do - - case(psb_dupl_add_) - ! Add - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= a%m)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nnz,a%ia1) - i2 = i1 - do - if (i2+1 > nnz) exit - if (a%ia1(i2+1) /= a%ia1(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (a%ia1(i1-1) /= a%ia1(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,a%ia2(i1:i2)) - if (ip>0) then - a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) - else - info = i - return - end if - end if - end do - - case default - info = -3 - if (debug_level >= psb_debug_serial_) & - & write(debug_unit,*) trim(name),& - & ': Duplicate handling: ',dupl - end select - - end if - - end subroutine d_coo_srch_upd - - subroutine d_jad_srch_upd(nz,ia,ja,val,nza,a,imin,imax,jmin,jmax,nzl,info,gtl,ng) - - use psb_spmat_type - use psb_const_mod - use psb_realloc_mod - use psb_string_mod - use psb_serial_mod - - implicit none - - type(psb_dspmat_type), intent(inout), target :: a - integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl - integer, intent(in) :: ia(:),ja(:) - integer, intent(inout) :: nza - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: ng,gtl(:) - - integer, pointer :: ia1(:), ia2(:), ia3(:),& - & ja_(:), ka_(:) - integer, allocatable :: indices(:), blks(:), rows(:) - integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ngr,& - & i,j,nr,dupl, ii, ir, ic - - info = 0 - dupl = psb_sp_getifld(psb_dupl_,a,info) - - png = a%ia2(1) ! points to the number of blocks - pia = a%ia2(2) ! points to the beginning of ia(3,png) - pja = a%ia2(3) ! points to the beginning of ja(:) - - ngr = a%ia2(png) ! the number of blocks - ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk - ka_ => a%ia1(:) ! the array containing the column indices - ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index - ! of each block - ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. - ! in ja to the first jad column - ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. - ! in ja to the first csr column - - - if (a%pl(1) /= 0) then - - if (present(gtl)) then - if (.not.present(ng)) then - info = -1 - return - endif - - - allocate(rows(nz),stat=info) - if (info /= 0) then - info = -4010 - return - endif - j=0 - do i=1,nz - ir = ia(i) - if ((ir >=1).and.(ir<=ng)) then - ir = gtl(ir) - j = j + 1 - rows(j) = ir - endif - enddo - call psb_msort_unique(rows(1:j),nr) - allocate(indices(nr),blks(nr),stat=info) - if (info /= 0) then - info = -4010 - return - endif - - do i=1,nr - indices(i)=a%pl(rows(i)) - j=0 - blkfnd_gtl: do - j=j+1 - if(ia1(j) == indices(i)) then - blks(i)=j - ipx = ia1(j) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j)+rb - exit blkfnd_gtl - else if(ia1(j) > indices(i)) then - blks(i)=j-1 - ipx = ia1(j-1) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j-1)+rb - exit blkfnd_gtl - end if - end do blkfnd_gtl - end do - - - ! cycle over elements - update_gtl: do ii=1,nz - ir = ia(ii) - ic = ja(ii) - if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then - ir = gtl(ir) - if ((ir > 0).and.(ir <= a%m)) then - ic = gtl(ic) - i = psb_ibsrch(ir,nr,rows) - - ! find which block the row belongs to - blk = blks(i) - - ! extract first part of the row from the jad block - ipx = ia1(blk) ! the first row index of the block - k_pt= ia2(blk) ! the pointer to the beginning of a column in ja - rb = indices(i)-ipx ! the row offset within the block - npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block - - do col = ia2(blk), ia3(blk)-1 - if (ic == ka_(ja_(col)+rb)) then - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - a%aspk(ja_(col)+rb) = val(ii) - case(psb_dupl_add_) - a%aspk(ja_(col)+rb) = a%aspk(ja_(col)+rb) + val(ii) - end select - cycle update_gtl - endif - end do - - ! extract second part of the row from the csr tail just in case - row=ia3(blk)+rb - do j=ja_(row), ja_(row+1)-1 - if (ic == ka_(j)) then - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - a%aspk(j) = val(ii) - case(psb_dupl_add_) - a%aspk(j) = a%aspk(j) + val(ii) - end select - cycle update_gtl - endif - end do - end if - end if - end do update_gtl - - else - - allocate(rows(nz),stat=info) - if (info /= 0) then - info = -4010 - return - endif - j=0 - do i=1,nz - ir = ia(i) - if ((ir >=1).and.(ir<=a%m)) then - j = j + 1 - rows(j) = ir - endif - enddo - call psb_msort_unique(rows(1:j),nr) - allocate(indices(nr),blks(nr),stat=info) - if (info /= 0) then - info = -4010 - return - endif - - do i=1,nr - indices(i)=a%pl(rows(i)) - j=0 - blkfnd: do - j=j+1 - if(ia1(j) == indices(i)) then - blks(i)=j - ipx = ia1(j) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j)+rb - exit blkfnd - else if(ia1(j) > indices(i)) then - blks(i)=j-1 - ipx = ia1(j-1) ! the first row index of the block - rb = indices(i)-ipx ! the row offset within the block - row = ia3(j-1)+rb - exit blkfnd - end if - end do blkfnd - end do - - - ! cycle over elements - update: do ii=1,nz - ir = ia(ii) - ic = ja(ii) - if ((ir >=1).and.(ir<=a%m).and.(ic>=1).and.(ic<=a%k)) then - i = psb_ibsrch(ir,nr,rows) - ! find which block the row belongs to - blk = blks(i) - - ! extract first part of the row from the jad block - ipx = ia1(blk) ! the first row index of the block - k_pt= ia2(blk) ! the pointer to the beginning of a column in ja - rb = indices(i)-ipx ! the row offset within the block - npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block - - do col = ia2(blk), ia3(blk)-1 - if (ic == ka_(ja_(col)+rb)) then - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - a%aspk(ja_(col)+rb) = val(ii) - case(psb_dupl_add_) - a%aspk(ja_(col)+rb) = a%aspk(ja_(col)+rb) + val(ii) - end select - cycle update - endif - end do - - ! extract second part of the row from the csr tail just in case - row=ia3(blk)+rb - do j=ja_(row), ja_(row+1)-1 - if (ic == ka_(j)) then - select case(dupl) - case(psb_dupl_ovwrt_,psb_dupl_err_) - a%aspk(j) = val(ii) - case(psb_dupl_add_) - a%aspk(j) = a%aspk(j) + val(ii) - end select - cycle update - endif - end do - end if - end do update - - end if - else - ! There might be some problems - info=134 - end if - - end subroutine d_jad_srch_upd - +!!$ subroutine d_csr_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_realloc_mod +!!$ use psb_string_mod +!!$ use psb_serial_mod +!!$ use psb_error_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ integer, intent(inout) :: nza +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: ng,gtl(:) +!!$ +!!$ integer :: debug_level, debug_unit +!!$ character(len=20) :: name='d_csr_srch_upd' +!!$ integer :: i,ir,ic, ilr, ilc, ip, & +!!$ & i1,i2,nc,lb,ub,m,dupl +!!$ +!!$ info = 0 +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ dupl = psb_sp_getifld(psb_dupl_,a,info) +!!$ +!!$ if (present(gtl)) then +!!$ if (.not.present(ng)) then +!!$ info = -1 +!!$ return +!!$ endif +!!$ +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ +!!$ ilr = -1 +!!$ ilc = -1 +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ ic = gtl(ic) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ i1 = a%ia2(ir) +!!$ i2 = a%ia2(ir+1) +!!$ nc=i2-i1 +!!$ +!!$ ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = val(i) +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Was searching ',ic,' in: ',i1,i2,& +!!$ & ' : ',a%ia1(i1:i2-1) +!!$ info = i +!!$ return +!!$ end if +!!$ +!!$ else +!!$ +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case(psb_dupl_add_) +!!$ ! Add +!!$ ilr = -1 +!!$ ilc = -1 +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ ic = gtl(ic) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ i1 = a%ia2(ir) +!!$ i2 = a%ia2(ir+1) +!!$ nc = i2-i1 +!!$ ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Was searching ',ic,' in: ',i1,i2,& +!!$ & ' : ',a%ia1(i1:i2-1) +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ +!!$ else +!!$ +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ +!!$ ilr = -1 +!!$ ilc = -1 +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ i1 = a%ia2(ir) +!!$ i2 = a%ia2(ir+1) +!!$ nc=i2-i1 +!!$ +!!$ ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = val(i) +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Was searching ',ic,' in: ',i1,i2,& +!!$ & ' : ',a%ia1(i1:i2-1) +!!$ info = i +!!$ return +!!$ end if +!!$ +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ +!!$ end do +!!$ +!!$ case(psb_dupl_add_) +!!$ ! Add +!!$ ilr = -1 +!!$ ilc = -1 +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ i1 = a%ia2(ir) +!!$ i2 = a%ia2(ir+1) +!!$ nc = i2-i1 +!!$ ip = psb_ibsrch(ic,nc,a%ia1(i1:i2-1)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ end if +!!$ +!!$ end subroutine d_csr_srch_upd +!!$ +!!$ subroutine d_coo_srch_upd(nz,ia,ja,val,nza,a,& +!!$ & imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_realloc_mod +!!$ use psb_string_mod +!!$ use psb_serial_mod +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout) :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ integer, intent(inout) :: nza +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: ng,gtl(:) +!!$ integer :: i,ir,ic, ilr, ilc, ip, & +!!$ & i1,i2,nc,nnz,dupl +!!$ integer :: debug_level, debug_unit +!!$ character(len=20) :: name='d_coo_srch_upd' +!!$ +!!$ info = 0 +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ dupl = psb_sp_getifld(psb_dupl_,a,info) +!!$ +!!$ if (psb_sp_getifld(psb_srtd_,a,info) /= psb_isrtdcoo_) then +!!$ info = -4 +!!$ return +!!$ end if +!!$ +!!$ ilr = -1 +!!$ ilc = -1 +!!$ nnz = psb_sp_getifld(psb_nnz_,a,info) +!!$ +!!$ +!!$ if (present(gtl)) then +!!$ if (.not.present(ng)) then +!!$ info = -1 +!!$ return +!!$ endif +!!$ +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ ic = gtl(ic) +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia1) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia1(i2+1) /= a%ia1(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia1(i1-1) /= a%ia1(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ia2(i1:i2)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ endif +!!$ end if +!!$ end do +!!$ case(psb_dupl_add_) +!!$ ! Add +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ ic = gtl(ic) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia1) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia1(i2+1) /= a%ia1(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia1(i1-1) /= a%ia1(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ia2(i1:i2)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ else +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Discarding row that does not belong to us.' +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ +!!$ else +!!$ +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia1) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia1(i2+1) /= a%ia1(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia1(i1-1) /= a%ia1(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ia2(i1:i2)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case(psb_dupl_add_) +!!$ ! Add +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia1) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia1(i2+1) /= a%ia1(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia1(i1-1) /= a%ia1(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ia2(i1:i2)) +!!$ if (ip>0) then +!!$ a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case default +!!$ info = -3 +!!$ if (debug_level >= psb_debug_serial_) & +!!$ & write(debug_unit,*) trim(name),& +!!$ & ': Duplicate handling: ',dupl +!!$ end select +!!$ +!!$ end if +!!$ +!!$ end subroutine d_coo_srch_upd +!!$ +!!$ subroutine d_jad_srch_upd(nz,ia,ja,val,nza,a,imin,imax,jmin,jmax,nzl,info,gtl,ng) +!!$ +!!$ use psb_spmat_type +!!$ use psb_const_mod +!!$ use psb_realloc_mod +!!$ use psb_string_mod +!!$ use psb_serial_mod +!!$ +!!$ implicit none +!!$ +!!$ type(psb_dspmat_type), intent(inout), target :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ integer, intent(inout) :: nza +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: ng,gtl(:) +!!$ +!!$ integer, pointer :: ia1(:), ia2(:), ia3(:),& +!!$ & ja_(:), ka_(:) +!!$ integer, allocatable :: indices(:), blks(:), rows(:) +!!$ integer :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ngr,& +!!$ & i,j,nr,dupl, ii, ir, ic +!!$ +!!$ info = 0 +!!$ dupl = psb_sp_getifld(psb_dupl_,a,info) +!!$ +!!$ png = a%ia2(1) ! points to the number of blocks +!!$ pia = a%ia2(2) ! points to the beginning of ia(3,png) +!!$ pja = a%ia2(3) ! points to the beginning of ja(:) +!!$ +!!$ ngr = a%ia2(png) ! the number of blocks +!!$ ja_ => a%ia2(pja:) ! the array containing the pointers to ka and aspk +!!$ ka_ => a%ia1(:) ! the array containing the column indices +!!$ ia1 => a%ia2(pia:pja-1:3) ! the array containing the first row index +!!$ ! of each block +!!$ ia2 => a%ia2(pia+1:pja-1:3) ! the array containing a pointer to the pos. +!!$ ! in ja to the first jad column +!!$ ia3 => a%ia2(pia+2:pja-1:3) ! the array containing a pointer to the pos. +!!$ ! in ja to the first csr column +!!$ +!!$ +!!$ if (a%pl(1) /= 0) then +!!$ +!!$ if (present(gtl)) then +!!$ if (.not.present(ng)) then +!!$ info = -1 +!!$ return +!!$ endif +!!$ +!!$ +!!$ allocate(rows(nz),stat=info) +!!$ if (info /= 0) then +!!$ info = -4010 +!!$ return +!!$ endif +!!$ j=0 +!!$ do i=1,nz +!!$ ir = ia(i) +!!$ if ((ir >=1).and.(ir<=ng)) then +!!$ ir = gtl(ir) +!!$ j = j + 1 +!!$ rows(j) = ir +!!$ endif +!!$ enddo +!!$ call psb_msort_unique(rows(1:j),nr) +!!$ allocate(indices(nr),blks(nr),stat=info) +!!$ if (info /= 0) then +!!$ info = -4010 +!!$ return +!!$ endif +!!$ +!!$ do i=1,nr +!!$ indices(i)=a%pl(rows(i)) +!!$ j=0 +!!$ blkfnd_gtl: do +!!$ j=j+1 +!!$ if(ia1(j) == indices(i)) then +!!$ blks(i)=j +!!$ ipx = ia1(j) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j)+rb +!!$ exit blkfnd_gtl +!!$ else if(ia1(j) > indices(i)) then +!!$ blks(i)=j-1 +!!$ ipx = ia1(j-1) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j-1)+rb +!!$ exit blkfnd_gtl +!!$ end if +!!$ end do blkfnd_gtl +!!$ end do +!!$ +!!$ +!!$ ! cycle over elements +!!$ update_gtl: do ii=1,nz +!!$ ir = ia(ii) +!!$ ic = ja(ii) +!!$ if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then +!!$ ir = gtl(ir) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ ic = gtl(ic) +!!$ i = psb_ibsrch(ir,nr,rows) +!!$ +!!$ ! find which block the row belongs to +!!$ blk = blks(i) +!!$ +!!$ ! extract first part of the row from the jad block +!!$ ipx = ia1(blk) ! the first row index of the block +!!$ k_pt= ia2(blk) ! the pointer to the beginning of a column in ja +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block +!!$ +!!$ do col = ia2(blk), ia3(blk)-1 +!!$ if (ic == ka_(ja_(col)+rb)) then +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ a%aspk(ja_(col)+rb) = val(ii) +!!$ case(psb_dupl_add_) +!!$ a%aspk(ja_(col)+rb) = a%aspk(ja_(col)+rb) + val(ii) +!!$ end select +!!$ cycle update_gtl +!!$ endif +!!$ end do +!!$ +!!$ ! extract second part of the row from the csr tail just in case +!!$ row=ia3(blk)+rb +!!$ do j=ja_(row), ja_(row+1)-1 +!!$ if (ic == ka_(j)) then +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ a%aspk(j) = val(ii) +!!$ case(psb_dupl_add_) +!!$ a%aspk(j) = a%aspk(j) + val(ii) +!!$ end select +!!$ cycle update_gtl +!!$ endif +!!$ end do +!!$ end if +!!$ end if +!!$ end do update_gtl +!!$ +!!$ else +!!$ +!!$ allocate(rows(nz),stat=info) +!!$ if (info /= 0) then +!!$ info = -4010 +!!$ return +!!$ endif +!!$ j=0 +!!$ do i=1,nz +!!$ ir = ia(i) +!!$ if ((ir >=1).and.(ir<=a%m)) then +!!$ j = j + 1 +!!$ rows(j) = ir +!!$ endif +!!$ enddo +!!$ call psb_msort_unique(rows(1:j),nr) +!!$ allocate(indices(nr),blks(nr),stat=info) +!!$ if (info /= 0) then +!!$ info = -4010 +!!$ return +!!$ endif +!!$ +!!$ do i=1,nr +!!$ indices(i)=a%pl(rows(i)) +!!$ j=0 +!!$ blkfnd: do +!!$ j=j+1 +!!$ if(ia1(j) == indices(i)) then +!!$ blks(i)=j +!!$ ipx = ia1(j) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j)+rb +!!$ exit blkfnd +!!$ else if(ia1(j) > indices(i)) then +!!$ blks(i)=j-1 +!!$ ipx = ia1(j-1) ! the first row index of the block +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ row = ia3(j-1)+rb +!!$ exit blkfnd +!!$ end if +!!$ end do blkfnd +!!$ end do +!!$ +!!$ +!!$ ! cycle over elements +!!$ update: do ii=1,nz +!!$ ir = ia(ii) +!!$ ic = ja(ii) +!!$ if ((ir >=1).and.(ir<=a%m).and.(ic>=1).and.(ic<=a%k)) then +!!$ i = psb_ibsrch(ir,nr,rows) +!!$ ! find which block the row belongs to +!!$ blk = blks(i) +!!$ +!!$ ! extract first part of the row from the jad block +!!$ ipx = ia1(blk) ! the first row index of the block +!!$ k_pt= ia2(blk) ! the pointer to the beginning of a column in ja +!!$ rb = indices(i)-ipx ! the row offset within the block +!!$ npg = ja_(k_pt+1)-ja_(k_pt) ! the number of rows in the block +!!$ +!!$ do col = ia2(blk), ia3(blk)-1 +!!$ if (ic == ka_(ja_(col)+rb)) then +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ a%aspk(ja_(col)+rb) = val(ii) +!!$ case(psb_dupl_add_) +!!$ a%aspk(ja_(col)+rb) = a%aspk(ja_(col)+rb) + val(ii) +!!$ end select +!!$ cycle update +!!$ endif +!!$ end do +!!$ +!!$ ! extract second part of the row from the csr tail just in case +!!$ row=ia3(blk)+rb +!!$ do j=ja_(row), ja_(row+1)-1 +!!$ if (ic == ka_(j)) then +!!$ select case(dupl) +!!$ case(psb_dupl_ovwrt_,psb_dupl_err_) +!!$ a%aspk(j) = val(ii) +!!$ case(psb_dupl_add_) +!!$ a%aspk(j) = a%aspk(j) + val(ii) +!!$ end select +!!$ cycle update +!!$ endif +!!$ end do +!!$ end if +!!$ end do update +!!$ +!!$ end if +!!$ else +!!$ ! There might be some problems +!!$ info=134 +!!$ end if +!!$ +!!$ end subroutine d_jad_srch_upd +!!$ subroutine c_csr_srch_upd(nz,ia,ja,val,nza,a,& & imin,imax,jmin,jmax,nzl,info,gtl,ng) diff --git a/base/tools/psb_dcdbldext.F90 b/base/tools/psb_dcdbldext.F90 index 12d8a18a..539f2a9d 100644 --- a/base/tools/psb_dcdbldext.F90 +++ b/base/tools/psb_dcdbldext.F90 @@ -67,6 +67,7 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) use psb_error_mod use psb_penv_mod use psb_realloc_mod + use psb_mat_mod use psi_mod #ifdef MPI_MOD use mpi @@ -78,14 +79,14 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) ! .. Array Arguments .. integer, intent(in) :: novr - Type(psb_dspmat_type), Intent(in) :: a + Type(psb_d_sparse_mat), Intent(in) :: a Type(psb_desc_type), Intent(in), target :: desc_a Type(psb_desc_type), Intent(out) :: desc_ov integer, intent(out) :: info integer, intent(in),optional :: extype ! .. Local Scalars .. - Integer :: i, j, np, me,m,nnzero,& + Integer :: i, j, np, me,m,& & ictxt, lovr, lworks,lworkr, n_row,n_col, int_err(5),& & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo Integer :: counter,counter_h, counter_o, counter_e,idx,gidx,proc,n_elem_recv,& @@ -94,7 +95,7 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) & idxr, idxs, iszr, iszs, nxch, nsnd, nrcv,lidx, extype_ integer :: icomm, err_act - type(psb_dspmat_type) :: blk + integer, allocatable :: irow(:), icol(:) Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:), orig_ovr(:) Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& & t_halo_out(:),temp(:),maskr(:) @@ -122,7 +123,6 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) extype_ = psb_ovt_xhal_ endif m = psb_cd_get_local_rows(desc_a) - nnzero = Size(a%aspk) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) nhalo = n_col-m @@ -169,7 +169,7 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) ! LOVR= (NNZ/NROW)*N_HALO*NOVR This assumes that the local average ! nonzeros per row is the same as the global. ! - nztot = psb_sp_get_nnzeros(a) + nztot = a%get_nzeros() if (nztot>0) then lovr = ((nztot+m-1)/m)*nhalo*novr lworks = ((nztot+m-1)/m)*nhalo @@ -210,16 +210,6 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) goto 9999 end if - call psb_sp_all(blk,max(lworks,lworkr),info) - if (info /= 0) then - info=4010 - ch_err='psb_sp_all' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - blk%fida='COO' - Allocate(orig_ovr(l_tmp_ovr_idx),tmp_ovr_idx(l_tmp_ovr_idx),& & tmp_halo(l_tmp_halo), halo(size(desc_a%halo_index)),stat=info) if (info /= 0) then @@ -414,35 +404,20 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) ! Prepare to exchange the halo rows with the other proc. ! If (i_ovr <= (novr)) Then - n_elem = psb_sp_get_nnz_row(idx,a) - call psb_ensure_size((idxs+tot_elem+n_elem),works,info) + call a%csget(idx,idx,n_elem,irow,icol,info) if (info /= 0) then info=4010 - call psb_errpush(info,name,a_err='psb_ensure_size') + call psb_errpush(info,name,a_err='csget') goto 9999 end if - - If((n_elem) > size(blk%ia2)) Then - isz = max((3*size(blk%ia2))/2,(n_elem)) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,'Realloc blk',isz - call psb_sp_reall(blk,isz,info) - if (info /= 0) then - info=4010 - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - End If - - call psb_sp_getblk(idx,a,blk,info) + + call psb_ensure_size((idxs+tot_elem+n_elem),works,info) if (info /= 0) then info=4010 - ch_err='psb_sp_getblk' - call psb_errpush(info,name,a_err=ch_err) + call psb_errpush(info,name,a_err='psb_ensure_size') goto 9999 end if - call psb_map_l2g(blk%ia2(1:n_elem),& + call psb_map_l2g(icol(1:n_elem),& & works(idxs+tot_elem+1:idxs+tot_elem+n_elem),& & desc_ov%idxmap,info) tot_elem=tot_elem+n_elem @@ -734,15 +709,21 @@ Subroutine psb_dcdbldext(a,desc_a,novr,desc_ov,info, extype) end if call psb_icdasb(desc_ov,info,ext_hv=.true.) - - call psb_cd_set_ovl_asb(desc_ov,info) - - if (info == 0) call psb_sp_free(blk,info) if (info /= 0) then - ch_err='sp_free' - call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + call psb_errpush(4010,name,a_err='icdasdb') goto 9999 end if + + call psb_cd_set_ovl_asb(desc_ov,info) + + if (info == 0) then + if (allocated(irow)) deallocate(irow,stat=info) + if ((info ==0).and.allocated(icol)) deallocate(icol,stat=info) + if (info /= 0) then + call psb_errpush(4013,name,a_err='deallocate',i_err=(/info,0,0,0,0/)) + goto 9999 + end if + end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),': end' diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 5ff153f3..feb63589 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -44,8 +44,6 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) use psb_descriptor_type - use psb_spmat_type - use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index cf9c4a82..50e9d09a 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -50,8 +50,6 @@ ! subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) use psb_descriptor_type - use psb_spmat_type - use psb_serial_mod use psb_const_mod use psi_mod use psb_error_mod diff --git a/base/tools/psb_dspfree.f90 b/base/tools/psb_dspfree.f90 index 6e1147e7..3e5a8161 100644 --- a/base/tools/psb_dspfree.f90 +++ b/base/tools/psb_dspfree.f90 @@ -42,8 +42,6 @@ subroutine psb_dspfree(a, desc_a,info) !...free sparse matrix structure... use psb_descriptor_type - use psb_spmat_type - use psb_serial_mod use psb_const_mod use psb_error_mod use psb_d_mat_mod @@ -73,13 +71,6 @@ subroutine psb_dspfree(a, desc_a,info) !...deallocate a.... call a%free() - -!!$ if(info /= 0) then -!!$ info=2045 -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if - call psb_erractionrestore(err_act) return diff --git a/base/tools/psb_dsphalo.F90 b/base/tools/psb_dsphalo.F90 index b2fbd359..62282f9d 100644 --- a/base/tools/psb_dsphalo.F90 +++ b/base/tools/psb_dsphalo.F90 @@ -60,7 +60,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& & rowscale,colscale,outfmt,data) use psb_const_mod - use psb_serial_mod + use psb_string_mod + use psb_mat_mod use psb_descriptor_type use psb_realloc_mod use psb_tools_mod, psb_protect_name => psb_dsphalo @@ -74,8 +75,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& include 'mpif.h' #endif - Type(psb_dspmat_type),Intent(in) :: a - Type(psb_dspmat_type),Intent(inout) :: blk + Type(psb_d_sparse_mat),Intent(in) :: a + Type(psb_d_sparse_mat),Intent(inout) :: blk Type(psb_desc_type),Intent(in), target :: desc_a integer, intent(out) :: info logical, optional, intent(in) :: rowcnv,colcnv,rowscale,colscale @@ -90,6 +91,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Integer, allocatable :: sdid(:,:), brvindx(:),rvid(:,:), & & rvsz(:), bsdindx(:),sdsz(:), iasnd(:), jasnd(:) real(psb_dpk_), allocatable :: valsnd(:) + type(psb_d_coo_sparse_mat), allocatable :: acoo integer, pointer :: idxv(:) logical :: rowcnv_,colcnv_,rowscale_,colscale_ character(len=5) :: outfmt_ @@ -144,7 +146,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Call psb_info(ictxt, me, np) Allocate(sdid(np,3),rvid(np,3),brvindx(np+1),& - & rvsz(np),sdsz(np),bsdindx(np+1),stat=info) + & rvsz(np),sdsz(np),bsdindx(np+1), acoo,stat=info) if (info /= 0) then info=4000 @@ -181,8 +183,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& idx = 0 idxs = 0 idxr = 0 - blk%k = a%k - blk%m = 0 + call acoo%allocate(0,a%get_ncols(),info) ! For all rows in the halo descriptor, extract and send/receive. Do proc=idxv(counter) @@ -193,13 +194,11 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& tot_elem = 0 Do j=0,n_el_send-1 idx = idxv(counter+psb_elem_send_+j) - n_elem = psb_sp_get_nnz_row(idx,a) + n_elem = a%get_nz_row(idx) tot_elem = tot_elem+n_elem Enddo sdsz(proc+1) = tot_elem - - blk%m = blk%m + n_el_recv - + call acoo%set_nrows(acoo%get_nrows() + n_el_recv) counter = counter+n_el_send+3 Enddo @@ -229,9 +228,9 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Enddo iszr=sum(rvsz) - call psb_sp_all(blk,max(iszr,1),info) + call acoo%reallocate(max(iszr,1)) if (debug_level >= psb_debug_outer_)& - & write(debug_unit,*) me,' ',trim(name),': Sizes:',size(blk%ia1),size(blk%ia2),& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),& & ' Send:',sdsz(:),' Receive:',rvsz(:) if (info /= 0) then info=4010 @@ -260,9 +259,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& Do j=0,n_el_send-1 idx = idxv(counter+psb_elem_send_+j) - n_elem = psb_sp_get_nnz_row(idx,a) - - call psb_sp_getrow(idx,a,ngtz,iasnd,jasnd,valsnd,info,& + n_elem = a%get_nz_row(idx) + call a%csget(idx,idx,ngtz,iasnd,jasnd,valsnd,info,& & append=.true.,nzin=tot_elem) if (info /= 0) then info=4010 @@ -272,9 +270,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& end if tot_elem=tot_elem+n_elem Enddo - ipx = ipx + 1 - counter = counter+n_el_send+3 Enddo nz = tot_elem @@ -290,11 +286,11 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& call mpi_alltoallv(valsnd,sdsz,bsdindx,mpi_double_precision,& - & blk%aspk,rvsz,brvindx,mpi_double_precision,icomm,info) + & acoo%val,rvsz,brvindx,mpi_double_precision,icomm,info) call mpi_alltoallv(iasnd,sdsz,bsdindx,mpi_integer,& - & blk%ia1,rvsz,brvindx,mpi_integer,icomm,info) + & acoo%ia,rvsz,brvindx,mpi_integer,icomm,info) call mpi_alltoallv(jasnd,sdsz,bsdindx,mpi_integer,& - & blk%ia2,rvsz,brvindx,mpi_integer,icomm,info) + & acoo%ja,rvsz,brvindx,mpi_integer,icomm,info) if (info /= 0) then info=4010 ch_err='mpi_alltoallv' @@ -305,8 +301,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& ! ! Convert into local numbering ! - if (rowcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I') - if (colcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') + if (rowcnv_) call psb_glob_to_loc(acoo%ia(1:iszr),desc_a,info,iact='I') + if (colcnv_) call psb_glob_to_loc(acoo%ja(1:iszr),desc_a,info,iact='I') if (info /= 0) then info=4010 @@ -316,21 +312,21 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& end if l1 = 0 - blk%m=0 + call acoo%set_nrows(0) ! irmin = huge(irmin) icmin = huge(icmin) irmax = 0 icmax = 0 Do i=1,iszr - r=(blk%ia1(i)) - k=(blk%ia2(i)) + r=(acoo%ia(i)) + k=(acoo%ja(i)) ! Just in case some of the conversions were out-of-range If ((r>0).and.(k>0)) Then l1=l1+1 - blk%aspk(l1) = blk%aspk(i) - blk%ia1(l1) = r - blk%ia2(l1) = k + acoo%val(l1) = acoo%val(i) + acoo%ia(l1) = r + acoo%ja(l1) = k irmin = min(irmin,r) irmax = max(irmax,r) icmin = min(icmin,k) @@ -338,37 +334,28 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rowcnv,colcnv,& End If Enddo if (rowscale_) then - blk%m = max(irmax-irmin+1,0) - blk%ia1(1:l1) = blk%ia1(1:l1) - irmin + 1 - else - blk%m = irmax + call acoo%set_nrows(max(irmax-irmin+1,0)) + acoo%ia(1:l1) = acoo%ia(1:l1) - irmin + 1 + else + call acoo%set_nrows(irmax) end if if (colscale_) then - blk%k = max(icmax-icmin+1,0) - blk%ia2(1:l1) = blk%ia2(1:l1) - icmin + 1 + call acoo%set_ncols(max(icmax-icmin+1,0)) + acoo%ja(1:l1) = acoo%ja(1:l1) - icmin + 1 else - blk%k = icmax + call acoo%set_ncols(icmax) end if + call acoo%set_nzeros(l1) - blk%fida = 'COO' - blk%infoa(psb_nnz_) = l1 - call psb_ensure_size(1,blk%pl,info) - if (info == 0) call psb_ensure_size(1,blk%pr,info) - if (info /= 0) then - info=4010 - ch_err='psb_ensure_size' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - blk%pl = 0 - blk%pr = 0 if (debug_level >= psb_debug_outer_)& & write(debug_unit,*) me,' ',trim(name),& - & ': End data exchange',counter,l1,blk%m + & ': End data exchange',counter,l1 + + call move_alloc(acoo,blk%a) ! Do we expect any duplicates to appear???? - call psb_spcnv(blk,info,afmt=outfmt_,dupl=psb_dupl_add_) + call blk%cscnv(info,type=outfmt_,dupl=psb_dupl_add_) if (info /= 0) then info=4010 ch_err='psb_spcnv' diff --git a/base/tools/psb_dspins.f90 b/base/tools/psb_dspins.f90 index a4074d5a..8ccd54cd 100644 --- a/base/tools/psb_dspins.f90 +++ b/base/tools/psb_dspins.f90 @@ -52,8 +52,6 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) use psb_tools_mod, psb_protect_name => psb_dspins use psb_descriptor_type - use psb_spmat_type - use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod @@ -241,17 +239,16 @@ end subroutine psb_dspins subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) use psb_tools_mod, psb_protect_name => psb_dspins_2desc use psb_descriptor_type - use psb_spmat_type - use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod + use psb_d_mat_mod implicit none !....parameters... type(psb_desc_type), intent(in) :: desc_ar type(psb_desc_type), intent(inout) :: desc_ac - type(psb_dspmat_type), intent(inout) :: a + type(psb_d_sparse_mat), intent(inout) :: a integer, intent(in) :: nz,ia(:),ja(:) real(kind=psb_dpk_), intent(in) :: val(:) integer, intent(out) :: info @@ -307,7 +304,6 @@ subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) end if if (nz==0) return - spstate = a%infoa(psb_state_) if (psb_is_bld_desc(desc_ac)) then allocate(ila(nz),jla(nz),stat=info) @@ -331,7 +327,7 @@ subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) nrow = psb_cd_get_local_rows(desc_ar) ncol = psb_cd_get_local_cols(desc_ac) - call psb_coins(nz,ila,jla,val,a,1,nrow,1,ncol,info) + call a%csput(nz,ila,jla,val,1,nrow,1,ncol,info) if (info /= 0) then info=4010 ch_err='psb_coins' diff --git a/base/tools/psb_dsprn.f90 b/base/tools/psb_dsprn.f90 index 2f7b25d1..817bea6c 100644 --- a/base/tools/psb_dsprn.f90 +++ b/base/tools/psb_dsprn.f90 @@ -44,7 +44,7 @@ Subroutine psb_dsprn(a, desc_a,info,clear) use psb_descriptor_type - use psb_spmat_type + use psb_mat_mod use psb_serial_mod use psb_const_mod use psb_error_mod @@ -53,7 +53,7 @@ Subroutine psb_dsprn(a, desc_a,info,clear) !....Parameters... Type(psb_desc_type), intent(in) :: desc_a - Type(psb_dspmat_type), intent(inout) :: a + Type(psb_d_sparse_mat), intent(inout) :: a integer, intent(out) :: info logical, intent(in), optional :: clear @@ -87,13 +87,8 @@ Subroutine psb_dsprn(a, desc_a,info,clear) call psb_errpush(info,name) goto 9999 endif - if (present(clear)) then - clear_ = clear - else - clear_ = .true. - end if - call psb_sp_reinit(a,info,clear=clear_) + call a%reinit(clear=clear) if (info /= 0) goto 9999 if (debug_level >= psb_debug_outer_) & diff --git a/base/tools/psb_linmap.f90 b/base/tools/psb_linmap.f90 index cb5ea07d..d1400c8b 100644 --- a/base/tools/psb_linmap.f90 +++ b/base/tools/psb_linmap.f90 @@ -113,7 +113,7 @@ function psb_d_linmap(map_kind,desc_X, desc_Y, map_X2Y, map_Y2X,iaggr,naggr) res implicit none type(psb_dlinmap_type) :: this type(psb_desc_type), target :: desc_X, desc_Y - type(psb_dspmat_type), intent(in) :: map_X2Y, map_Y2X + type(psb_d_sparse_mat), intent(in) :: map_X2Y, map_Y2X integer, intent(in) :: map_kind integer, intent(in), optional :: iaggr(:), naggr(:) ! @@ -171,8 +171,8 @@ function psb_d_linmap(map_kind,desc_X, desc_Y, map_X2Y, map_Y2X,iaggr,naggr) res info = 1 end select - if (info == 0) call psb_sp_clone(map_X2Y,this%map_X2Y,info) - if (info == 0) call psb_sp_clone(map_Y2X,this%map_Y2X,info) + if (info == 0) call psb_clone(map_X2Y,this%map_X2Y,info) + if (info == 0) call psb_clone(map_Y2X,this%map_Y2X,info) if (info == 0) call psb_realloc(psb_itd_data_size_,this%itd_data,info) if (info == 0) then call psb_set_map_kind(map_kind, this) @@ -182,8 +182,8 @@ function psb_d_linmap(map_kind,desc_X, desc_Y, map_X2Y, map_Y2X,iaggr,naggr) res return end if if (debug) then - write(0,*) trim(name),' forward map:',allocated(this%map_X2Y%aspk) - write(0,*) trim(name),' backward map:',allocated(this%map_Y2X%aspk) +!!$ write(0,*) trim(name),' forward map:',allocated(this%map_X2Y%aspk) +!!$ write(0,*) trim(name),' backward map:',allocated(this%map_Y2X%aspk) end if end function psb_d_linmap diff --git a/krylov/psb_dbicg.f90 b/krylov/psb_dbicg.f90 index 42e65ac7..b3bde503 100644 --- a/krylov/psb_dbicg.f90 +++ b/krylov/psb_dbicg.f90 @@ -98,12 +98,9 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_dbicg - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ parameters -!!$ type(psb_dspmat_type), intent(in) :: a type(psb_dprec_type), intent(in) :: prec type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: b(:) diff --git a/krylov/psb_dcg.F90 b/krylov/psb_dcg.F90 index 62e6ae21..4caee811 100644 --- a/krylov/psb_dcg.F90 +++ b/krylov/psb_dcg.F90 @@ -99,14 +99,11 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop,cond) use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_dcg - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ Parameters -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_dprec_type), Intent(in) :: prec Type(psb_desc_type), Intent(in) :: desc_a Real(psb_dpk_), Intent(in) :: b(:) diff --git a/krylov/psb_dcgs.f90 b/krylov/psb_dcgs.f90 index 92cc8dc3..831ed963 100644 --- a/krylov/psb_dcgs.f90 +++ b/krylov/psb_dcgs.f90 @@ -98,13 +98,10 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_dcgs - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ parameters -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a Type(psb_dprec_type), Intent(in) :: prec Real(psb_dpk_), Intent(in) :: b(:) diff --git a/krylov/psb_dcgstab.F90 b/krylov/psb_dcgstab.F90 index a3e9ca20..b211e513 100644 --- a/krylov/psb_dcgstab.F90 +++ b/krylov/psb_dcgstab.F90 @@ -98,13 +98,10 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop) use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_dcgstab - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ parameters -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_dprec_type), Intent(in) :: prec Type(psb_desc_type), Intent(in) :: desc_a Real(psb_dpk_), Intent(in) :: b(:) diff --git a/krylov/psb_dcgstabl.f90 b/krylov/psb_dcgstabl.f90 index 2076a80e..f15aacf2 100644 --- a/krylov/psb_dcgstabl.f90 +++ b/krylov/psb_dcgstabl.f90 @@ -107,13 +107,10 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_dcgstabl - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ parameters -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_dprec_type), Intent(in) :: prec Type(psb_desc_type), Intent(in) :: desc_a Real(psb_dpk_), Intent(in) :: b(:) diff --git a/krylov/psb_drgmres.f90 b/krylov/psb_drgmres.f90 index f9ae316f..5edcdbad 100644 --- a/krylov/psb_drgmres.f90 +++ b/krylov/psb_drgmres.f90 @@ -110,13 +110,10 @@ subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist use psb_base_mod use psb_prec_mod use psb_krylov_mod, psb_protect_name => psb_drgmres - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a -!!$ Parameters -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_dprec_type), Intent(in) :: prec Type(psb_desc_type), Intent(in) :: desc_a Real(psb_dpk_), Intent(in) :: b(:) diff --git a/krylov/psb_krylov_mod.f90 b/krylov/psb_krylov_mod.f90 index d6aab334..392818a7 100644 --- a/krylov/psb_krylov_mod.f90 +++ b/krylov/psb_krylov_mod.f90 @@ -60,9 +60,8 @@ Module psb_krylov_mod end subroutine psb_scg subroutine psb_dcg(a,prec,b,x,eps,& & desc_a,info,itmax,iter,err,itrace,istop,cond) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: b(:) @@ -124,9 +123,8 @@ Module psb_krylov_mod end subroutine psb_sbicg subroutine psb_dbicg(a,prec,b,x,eps,& & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: b(:) @@ -188,9 +186,8 @@ Module psb_krylov_mod end subroutine psb_scgstab subroutine psb_dcgstab(a,prec,b,x,eps,& & desc_a,info,itmax,iter,err,itrace,istop) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: b(:) @@ -252,9 +249,8 @@ Module psb_krylov_mod end subroutine psb_scgstabl Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& &itmax,iter,err, itrace,irst,istop) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a type(psb_dprec_type), intent(in) :: prec @@ -316,9 +312,8 @@ Module psb_krylov_mod end subroutine psb_srgmres Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace,irst,istop) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a type(psb_dprec_type), intent(in) :: prec @@ -380,9 +375,8 @@ Module psb_krylov_mod end subroutine psb_scgs subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& &itmax,iter,err,itrace,istop) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_mod, only : psb_dprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a type(psb_dprec_type), intent(in) :: prec @@ -585,7 +579,7 @@ contains ! BICGSTABL ! RGMRES ! - ! a - type(psb_dspmat_type) Input: sparse matrix containing A. + ! a - type(psb_d_sparse_mat) Input: sparse matrix containing A. ! prec - type(psb_dprec_type) Input: preconditioner ! b - real,dimension(:) Input: vector containing the ! right hand side B @@ -619,12 +613,10 @@ contains use psb_base_mod use psb_prec_mod,only : psb_sprec_type, psb_dprec_type, psb_cprec_type, psb_zprec_type - use psb_d_mat_mod type(psb_d_sparse_mat), intent(in) :: a character(len=*) :: method -!!$ Type(psb_dspmat_type), Intent(in) :: a Type(psb_desc_type), Intent(in) :: desc_a type(psb_dprec_type), intent(in) :: prec Real(psb_dpk_), Intent(in) :: b(:) @@ -1068,12 +1060,10 @@ contains subroutine psb_d_init_conv(methdname,stopc,trace,itmax,a,b,eps,desc_a,stopdat,info) use psb_base_mod - use psb_d_mat_mod implicit none type(psb_d_sparse_mat), intent(in) :: a character(len=*), intent(in) :: methdname integer, intent(in) :: stopc, trace,itmax -!!$ type(psb_dspmat_type), intent(in) :: a real(psb_dpk_), intent(in) :: b(:), eps type(psb_desc_type), intent(in) :: desc_a type(psb_itconv_type) :: stopdat diff --git a/prec/psb_dbjac_aply.f90 b/prec/psb_dbjac_aply.f90 index be65522c..23877eed 100644 --- a/prec/psb_dbjac_aply.f90 +++ b/prec/psb_dbjac_aply.f90 @@ -38,7 +38,6 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_d_mat_mod use psb_prec_mod, psb_protect_name => psb_dbjac_aply implicit none diff --git a/prec/psb_dbjac_bld.f90 b/prec/psb_dbjac_bld.f90 index 1006bb57..57a05ec9 100644 --- a/prec/psb_dbjac_bld.f90 +++ b/prec/psb_dbjac_bld.f90 @@ -30,7 +30,6 @@ !!$ !!$ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) - use psb_d_mat_mod use psb_base_mod use psb_prec_mod, psb_protect_name => psb_dbjac_bld implicit none diff --git a/prec/psb_ddiagsc_bld.f90 b/prec/psb_ddiagsc_bld.f90 index fedba61a..327b8770 100644 --- a/prec/psb_ddiagsc_bld.f90 +++ b/prec/psb_ddiagsc_bld.f90 @@ -32,7 +32,6 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) use psb_base_mod - use psb_d_mat_mod use psb_prec_mod, psb_protect_name => psb_ddiagsc_bld Implicit None diff --git a/prec/psb_dilu_fct.f90 b/prec/psb_dilu_fct.f90 index efb2365e..6aca9652 100644 --- a/prec/psb_dilu_fct.f90 +++ b/prec/psb_dilu_fct.f90 @@ -37,7 +37,6 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck) ! ! use psb_base_mod - use psb_d_mat_mod implicit none ! .. Scalar Arguments .. integer, intent(out) :: info diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index 05cd6d17..46da5e51 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -32,7 +32,6 @@ subroutine psb_dprecbld(aa,desc_a,p,info,upd) use psb_base_mod - use psb_d_mat_mod use psb_prec_mod, psb_protect_name => psb_dprecbld Implicit None diff --git a/prec/psb_prec_mod.f90 b/prec/psb_prec_mod.f90 index 5ca2f840..6639a785 100644 --- a/prec/psb_prec_mod.f90 +++ b/prec/psb_prec_mod.f90 @@ -45,8 +45,7 @@ module psb_prec_mod character, intent(in),optional :: upd end subroutine psb_sprecbld subroutine psb_dprecbld(a,desc_a,prec,info,upd) - use psb_d_mat_mod - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_type, only : psb_dprec_type implicit none type(psb_d_sparse_mat), intent(in), target :: a @@ -329,14 +328,12 @@ module psb_prec_mod real(psb_spk_), intent(inout) :: d(:) end subroutine psb_silu_fct subroutine psb_dilu_fct(a,l,u,d,info,blck) - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ - use psb_d_mat_mod + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat,& + & psb_d_csr_sparse_mat, psb_dpk_ integer, intent(out) :: info type(psb_d_sparse_mat),intent(in) :: a type(psb_d_csr_sparse_mat),intent(inout) :: l,u -!!$ type(psb_dspmat_type),intent(in) :: a -!!$ type(psb_dspmat_type),intent(inout) :: l,u - type(psb_dspmat_type),intent(in), optional, target :: blck + type(psb_d_sparse_mat),intent(in), optional, target :: blck real(psb_dpk_), intent(inout) :: d(:) end subroutine psb_dilu_fct subroutine psb_cilu_fct(a,l,u,d,info,blck) @@ -368,8 +365,7 @@ module psb_prec_mod character, intent(in) :: upd end subroutine psb_sbjac_bld subroutine psb_dbjac_bld(a,desc_a,p,upd,info) - use psb_d_mat_mod - use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + use psb_base_mod, only : psb_desc_type, psb_d_sparse_mat, psb_dpk_ use psb_prec_type, only : psb_dprec_type integer, intent(out) :: info type(psb_d_sparse_mat), intent(in), target :: a diff --git a/test/fileread/df_sample.f90 b/test/fileread/df_sample.f90 index bba9629f..4328cfdd 100644 --- a/test/fileread/df_sample.f90 +++ b/test/fileread/df_sample.f90 @@ -41,7 +41,7 @@ program df_sample character(len=40) :: kmethd, ptype, mtrx_file, rhs_file ! sparse matrices - type(psb_dspmat_type) :: a, aux_a + type(psb_d_sparse_mat) :: a, aux_a ! preconditioner data type(psb_dprec_type) :: prec @@ -129,7 +129,7 @@ program df_sample call psb_abort(ictxt) end if - m_problem = aux_a%m + m_problem = aux_a%get_nrows() call psb_bcast(ictxt,m_problem) ! At this point aux_b may still be unallocated @@ -182,7 +182,13 @@ program df_sample write(*,'("Partition type: graph")') write(*,'(" ")') ! write(0,'("Build type: graph")') - call build_mtpart(aux_a%m,aux_a%fida,aux_a%ia1,aux_a%ia2,np) + select type (aa=>a%a) + type is (psb_d_csr_sparse_mat) + call build_mtpart(aa%get_nrows(),aa%get_fmt(),aa%ja,aa%irp,np) + class default + write(0,*) 'Should never get here!' + call psb_abort(ictxt) + end select endif call psb_barrier(ictxt) call distr_mtpart(psb_root_,ictxt) diff --git a/test/fileread/runs/dfs.inp b/test/fileread/runs/dfs.inp index 1be6c6de..8027620d 100644 --- a/test/fileread/runs/dfs.inp +++ b/test/fileread/runs/dfs.inp @@ -1,11 +1,11 @@ 11 Number of inputs -sherman3.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or -sherman3_rhs1.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html +thm200x120.mtx sherman3.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or +NONE sherman3_rhs1.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html MM File format: MM: Matrix Market HB: Harwell-Boeing. BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG BJAC Preconditioner NONE DIAG BJAC CSR Storage format CSR COO JAD -0 IPART: Partition method 0: BLK 2: graph (with Metis) +2 IPART: Partition method 0: BLK 2: graph (with Metis) 2 ISTOPC 01000 ITMAX 01 ITRACE diff --git a/test/pargen/ppde.f90 b/test/pargen/ppde.f90 index 25bb26e9..873160ad 100644 --- a/test/pargen/ppde.f90 +++ b/test/pargen/ppde.f90 @@ -65,7 +65,6 @@ program ppde use psb_base_mod use psb_prec_mod use psb_krylov_mod - use psb_d_mat_mod implicit none ! input parameters diff --git a/util/psb_hbio_mod.f90 b/util/psb_hbio_mod.f90 index f1c243ae..d8977abd 100644 --- a/util/psb_hbio_mod.f90 +++ b/util/psb_hbio_mod.f90 @@ -328,7 +328,7 @@ contains subroutine dhb_read(a, iret, iunit, filename,b,g,x,mtitle) use psb_base_mod implicit none - type(psb_dspmat_type), intent(out) :: a + type(psb_d_sparse_mat), intent(out) :: a integer, intent(out) :: iret integer, optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename @@ -340,6 +340,8 @@ contains character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20 integer :: indcrd, ptrcrd, totcrd,& & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix + type(psb_d_csr_sparse_mat) :: acsr + type(psb_d_coo_sparse_mat) :: acoo integer :: ircode, i,nzr,infile, info character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)' character(len=*), parameter :: fmt11='(a3,11x,2i14)' @@ -369,27 +371,26 @@ contains read (infile,fmt=fmt10) mtitle_,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,& & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt if (rhscrd > 0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix - - call psb_sp_all(a,nnzero,nrow+1,nnzero,ircode) + + call acsr%allocate(nrow,ncol,nnzero) if (ircode /= 0 ) then write(0,*) 'Memory allocation failed' goto 993 end if + if (present(mtitle)) mtitle=mtitle_ - a%m = nrow - a%k = ncol - a%fida = 'CSR' - a%descra='G' - + if (psb_tolower(type(1:1)) == 'r') then if (psb_tolower(type(2:2)) == 'u') then - read (infile,fmt=ptrfmt) (a%ia2(i),i=1,nrow+1) - read (infile,fmt=indfmt) (a%ia1(i),i=1,nnzero) - if (valcrd > 0) read (infile,fmt=valfmt) (a%aspk(i),i=1,nnzero) - + read (infile,fmt=ptrfmt) (acsr%irp(i),i=1,nrow+1) + read (infile,fmt=indfmt) (acsr%ja(i),i=1,nnzero) + if (valcrd > 0) read (infile,fmt=valfmt) (acsr%val(i),i=1,nnzero) + + call a%mv_from(acsr) + if (present(b)) then if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then call psb_realloc(nrow,1,b,info) @@ -414,9 +415,9 @@ contains ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read - read (infile,fmt=ptrfmt) (a%ia2(i),i=1,nrow+1) - read (infile,fmt=indfmt) (a%ia1(i),i=1,nnzero) - if (valcrd > 0) read (infile,fmt=valfmt) (a%aspk(i),i=1,nnzero) + read (infile,fmt=ptrfmt) (acsr%irp(i),i=1,nrow+1) + read (infile,fmt=indfmt) (acsr%ja(i),i=1,nnzero) + if (valcrd > 0) read (infile,fmt=valfmt) (acsr%val(i),i=1,nnzero) if (present(b)) then @@ -438,23 +439,24 @@ contains endif endif - call psb_spcnv(a,ircode,afmt='csr') - if (ircode /= 0) goto 993 - - call psb_sp_reall(a,2*nnzero,ircode) + + call acoo%mv_from_fmt(acsr,info) + call acoo%reallocate(2*nnzero) ! A is now in COO format nzr = nnzero do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then + if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - a%aspk(nzr) = a%aspk(i) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) + acoo%val(nzr) = acoo%val(i) + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) end if end do - a%infoa(psb_nnz_) = nzr - call psb_spcnv(a,ircode,afmt='csr') - if (ircode /= 0) goto 993 + call acoo%set_nzeros(nzr) + call acoo%fix(ircode) + if (ircode==0) call a%mv_from(acoo) + if (ircode==0) call a%cscnv(ircode,type='csr') + if (ircode/=0) goto 993 else write(0,*) 'read_matrix: matrix type not yet supported' @@ -483,7 +485,7 @@ contains subroutine dhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle) use psb_base_mod implicit none - type(psb_dspmat_type), intent(in) :: a + type(psb_d_sparse_mat), intent(in) :: a integer, intent(out) :: iret character(len=*), optional, intent(in) :: mtitle integer, optional, intent(in) :: iunit @@ -540,11 +542,13 @@ contains key_ = 'PSBMAT00' endif - if (psb_toupper(a%fida) == 'CSR') then + + select type(aa=>a%a) + type is (psb_d_csr_sparse_mat) - nrow = a%m - ncol = a%k - nnzero = a%ia2(nrow+1)-1 + nrow = aa%get_nrows() + ncol = aa%get_ncols() + nnzero = aa%get_nzeros() neltvl = 0 @@ -583,19 +587,19 @@ contains write (iout,fmt=fmt10) mtitle_,key_,totcrd,ptrcrd,indcrd,valcrd,rhscrd,& & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt if (rhscrd > 0) write (iout,fmt=fmt11) rhstype,nrhs,nrhsix - write (iout,fmt=ptrfmt) (a%ia2(i),i=1,nrow+1) - write (iout,fmt=indfmt) (a%ia1(i),i=1,nnzero) - if (valcrd > 0) write (iout,fmt=valfmt) (a%aspk(i),i=1,nnzero) + write (iout,fmt=ptrfmt) (aa%irp(i),i=1,nrow+1) + write (iout,fmt=indfmt) (aa%ja(i),i=1,nnzero) + if (valcrd > 0) write (iout,fmt=valfmt) (aa%val(i),i=1,nnzero) if (rhscrd > 0) write (iout,fmt=rhsfmt) (rhs(i),i=1,nrow) if (present(g).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (g(i),i=1,nrow) if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow) - else + class default - write(0,*) 'format: ',a%fida,' not yet implemented' + write(0,*) 'format: ',a%get_fmt(),' not yet implemented' - endif + end select if (iout /= 6) close(iout) diff --git a/util/psb_mat_dist_mod.f90 b/util/psb_mat_dist_mod.f90 index 15718797..ca3b696c 100644 --- a/util/psb_mat_dist_mod.f90 +++ b/util/psb_mat_dist_mod.f90 @@ -551,11 +551,11 @@ contains implicit none ! parameters - type(psb_dspmat_type) :: a_glob + type(psb_d_sparse_mat) :: a_glob real(psb_dpk_) :: b_glob(:) integer :: ictxt type(psb_d_sparse_mat) :: a - real(psb_dpk_), allocatable :: b(:) + real(psb_dpk_), allocatable :: b(:) type(psb_desc_type) :: desc_a integer, intent(out) :: info integer, optional :: inroot @@ -599,22 +599,15 @@ contains end if call psb_info(ictxt, iam, np) if (iam == root) then - ! extract information from a_glob - if (a_glob%fida /= 'CSR') then - info=135 - ch_err='CSR' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - nrow = a_glob%m - ncol = a_glob%k + nrow = a_glob%get_nrows() + ncol = a_glob%get_ncols() if (nrow /= ncol) then write(0,*) 'a rectangular matrix ? ',nrow,ncol info=-1 call psb_errpush(info,name) goto 9999 endif - nnzero = size(a_glob%aspk) + nnzero = a_glob%get_nzeros() nrhs = 1 endif @@ -719,7 +712,7 @@ contains ll = 0 do i= i_count, j_count-1 - call psb_sp_getrow(i,a_glob,nz,& + call a_glob%csget(i,i,nz,& & irow,icol,val,info,nzin=ll,append=.true.) if (info /= 0) then if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then @@ -807,7 +800,7 @@ contains ll = 0 do i= i_count, i_count - call psb_sp_getrow(i,a_glob,nz,& + call a_glob%csget(i,i,nz,& & irow,icol,val,info,nzin=ll,append=.true.) if (info /= 0) then if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then diff --git a/util/psb_mmio_mod.f90 b/util/psb_mmio_mod.f90 index bd860c33..26339dd5 100644 --- a/util/psb_mmio_mod.f90 +++ b/util/psb_mmio_mod.f90 @@ -491,7 +491,7 @@ contains subroutine dmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none - type(psb_dspmat_type), intent(out) :: a + type(psb_d_sparse_mat), intent(out) :: a integer, intent(out) :: info integer, optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename @@ -499,6 +499,7 @@ contains character(1024) :: line integer :: nrow, ncol, nnzero integer :: ircode, i,nzr,infile + type(psb_d_coo_sparse_mat), allocatable :: acoo info = 0 @@ -534,46 +535,51 @@ contains if (line(1:1) /= '%') exit end do read(line,fmt=*) nrow,ncol,nnzero - + + allocate(acoo, stat=ircode) + if (ircode /= 0) goto 993 if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then - call psb_sp_all(nrow,ncol,a,nnzero,ircode) - a%fida = 'COO' - a%descra = 'G' - if (ircode /= 0) goto 993 + call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),a%aspk(i) + read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),acoo%val(i) end do - a%infoa(psb_nnz_) = nnzero - call psb_spcnv(a,ircode,afmt='csr') + call acoo%set_nzeros(nnzero) + call acoo%fix(info) + + call a%mv_from(acoo) + call a%cscnv(ircode,type='csr') else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read - call psb_sp_all(nrow,ncol,a,2*nnzero,ircode) - a%fida = 'COO' - a%descra = 'G' - if (ircode /= 0) goto 993 + call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),a%aspk(i) + read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),acoo%val(i) end do - nzr = nnzero do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then + if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - a%aspk(nzr) = a%aspk(i) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) + acoo%val(nzr) = acoo%val(i) + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) end if end do - a%infoa(psb_nnz_) = nzr - call psb_spcnv(a,ircode,afmt='csr') + call acoo%set_nzeros(nzr) + call acoo%fix(info) + + call a%mv_from(acoo) + call a%cscnv(ircode,type='csr') else write(0,*) 'read_matrix: matrix type not yet supported' info=904 end if + + if (infile/=5) close(infile) + + return ! open failed @@ -592,7 +598,7 @@ contains subroutine dmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none - type(psb_dspmat_type), intent(in) :: a + type(psb_d_sparse_mat), intent(in) :: a integer, intent(out) :: info character(len=*), intent(in) :: mtitle integer, optional, intent(in) :: iunit @@ -621,7 +627,7 @@ contains endif endif - call psb_csprt(iout,a,head=mtitle) + call a%print(iout,head=mtitle) if (iout /= 6) close(iout)