From 4ecc1b632d38f7cedd05c9f8359613b015a8f8a1 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 14 Sep 2009 13:55:43 +0000 Subject: [PATCH] psblas3: base/modules/psb_error_mod.F90 base/modules/psb_psblas_mod.f90 base/newserial/psbn_base_mat_mod.f03 base/newserial/psbn_d_base_mat_mod.f03 base/newserial/psbn_d_coo_impl.f03 base/newserial/psbn_d_csr_impl.f03 base/newserial/psbn_d_csr_mat_mod.f03 base/newserial/psbn_mat_mod.f03 base/psblas/psb_dnrmi.f90 base/psblas/psb_dspmm.f90 base/psblas/psb_dspsm.f90 base/tools/psb_dspalloc.f90 prec/psb_dbjac_aply.f90 prec/psb_dbjac_bld.f90 prec/psb_dilu_fct.f90 prec/psb_dprecbld.f90 prec/psb_dprecinit.f90 prec/psb_prec_mod.f90 prec/psb_prec_type.f90 test/pargen/ppde.f90 test/pargen/runs/ppde.inp Now both BJAC_BLD and CSSV work. Really! And initial performance is not too bad. Lots and lots of details to be fixed yet........... --- base/modules/psb_error_mod.F90 | 6 + base/modules/psb_psblas_mod.f90 | 6 +- base/newserial/psbn_base_mat_mod.f03 | 16 +- base/newserial/psbn_d_base_mat_mod.f03 | 356 ++++++++++++++++++++++++- base/newserial/psbn_d_coo_impl.f03 | 62 ++--- base/newserial/psbn_d_csr_impl.f03 | 247 +++++++++-------- base/newserial/psbn_d_csr_mat_mod.f03 | 18 +- base/newserial/psbn_mat_mod.f03 | 60 ++++- base/psblas/psb_dnrmi.f90 | 2 +- base/psblas/psb_dspmm.f90 | 10 +- base/psblas/psb_dspsm.f90 | 41 +-- base/tools/psb_dspalloc.f90 | 2 +- prec/psb_dbjac_aply.f90 | 1 + prec/psb_dbjac_bld.f90 | 68 +++-- prec/psb_dilu_fct.f90 | 109 +++----- prec/psb_dprecbld.f90 | 2 +- prec/psb_dprecinit.f90 | 12 +- prec/psb_prec_mod.f90 | 10 +- prec/psb_prec_type.f90 | 7 +- test/pargen/ppde.f90 | 7 +- test/pargen/runs/ppde.inp | 8 +- 21 files changed, 725 insertions(+), 325 deletions(-) diff --git a/base/modules/psb_error_mod.F90 b/base/modules/psb_error_mod.F90 index 4019065a..c45d68b4 100644 --- a/base/modules/psb_error_mod.F90 +++ b/base/modules/psb_error_mod.F90 @@ -386,9 +386,15 @@ contains case(30) write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1) write (0,'("current value is ",i0)')i_e_d(2) + case(31) + write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1) + write (0,'("current value is ",a)')a_e_d case(35) write (0,'("Size of input array argument n. ",i0," is invalid.")')i_e_d(1) write (0,'("Current value is ",i0)')i_e_d(2) + case(36) + write (0,'("Size of input array argument n. ",i0," must be ")')i_e_d(1) + write (0,'("at least ",i0)')i_e_d(2) case(40) write (0,'("input argument n. ",i0," has an invalid value")')i_e_d(1) write (0,'("current value is ",a)')a_e_d(2:2) diff --git a/base/modules/psb_psblas_mod.f90 b/base/modules/psb_psblas_mod.f90 index 8eaba543..80cb8c08 100644 --- a/base/modules/psb_psblas_mod.f90 +++ b/base/modules/psb_psblas_mod.f90 @@ -771,7 +771,8 @@ module psb_psblas_mod & diag, n, jx, jy, work) use psb_serial_mod use psb_descriptor_type - type(psb_dspmat_type), intent(in) :: t + use psbn_d_mat_mod + type(psbn_d_sparse_mat), intent(in) :: t real(psb_dpk_), intent(in) :: x(:,:) real(psb_dpk_), intent(inout) :: y(:,:) real(psb_dpk_), intent(in) :: alpha, beta @@ -787,7 +788,8 @@ module psb_psblas_mod & diag, work) use psb_serial_mod use psb_descriptor_type - type(psb_dspmat_type), intent(in) :: t + use psbn_d_mat_mod + type(psbn_d_sparse_mat), intent(in) :: t real(psb_dpk_), intent(in) :: x(:) real(psb_dpk_), intent(inout) :: y(:) real(psb_dpk_), intent(in) :: alpha, beta diff --git a/base/newserial/psbn_base_mat_mod.f03 b/base/newserial/psbn_base_mat_mod.f03 index 1aa3bc43..344bbde4 100644 --- a/base/newserial/psbn_base_mat_mod.f03 +++ b/base/newserial/psbn_base_mat_mod.f03 @@ -88,7 +88,8 @@ module psbn_base_mat_mod generic, public :: allocate => allocate_mnnz generic, public :: reallocate => reallocate_nz procedure, pass(a) :: print => sparse_print - + procedure, pass(a) :: sizeof + end type psbn_base_sparse_mat private :: set_nrows, set_ncols, set_dupl, set_state, & @@ -97,10 +98,21 @@ module psbn_base_mat_mod & get_nzeros, get_size, get_state, get_dupl, is_null, is_bld, & & is_upd, is_asb, is_sorted, is_upper, is_lower, is_triangle, & & is_unit, get_neigh, allocate_mn, allocate_mnnz, reallocate_nz, & - & free, sparse_print, get_fmt, trim + & free, sparse_print, get_fmt, trim, sizeof + contains + + + function sizeof(a) result(res) + implicit none + class(psbn_base_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + end function sizeof + + function get_fmt(a) result(res) implicit none class(psbn_base_sparse_mat), intent(in) :: a diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 39109db1..dc6cd054 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -9,13 +9,17 @@ module psbn_d_base_mat_mod generic, public :: csmm => d_base_csmm, d_base_csmv procedure, pass(a) :: d_base_cssv procedure, pass(a) :: d_base_cssm - generic, public :: cssm => d_base_cssm, d_base_cssv + generic, public :: base_cssm => d_base_cssm, d_base_cssv + procedure, pass(a) :: d_cssv + procedure, pass(a) :: d_cssm + generic, public :: cssm => d_cssm, d_cssv procedure, pass(a) :: d_scals procedure, pass(a) :: d_scal generic, public :: scal => d_scals, d_scal - procedure, pass(a) :: get_diag procedure, pass(a) :: csnmi + procedure, pass(a) :: get_diag procedure, pass(a) :: csput + procedure, pass(a) :: d_csgetrow procedure, pass(a) :: d_csgetblk generic, public :: csget => d_csgetrow, d_csgetblk @@ -34,7 +38,7 @@ module psbn_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 + & get_diag, csclip, d_cssv, d_cssm type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat @@ -73,6 +77,7 @@ module psbn_d_base_mat_mod procedure, pass(a) :: d_csgetrow => d_coo_csgetrow procedure, pass(a) :: print => d_coo_print procedure, pass(a) :: get_fmt => d_coo_get_fmt + procedure, pass(a) :: sizeof => d_coo_sizeof end type psbn_d_coo_sparse_mat @@ -82,7 +87,7 @@ module psbn_d_base_mat_mod & d_fix_coo, d_coo_free, d_coo_print, d_coo_get_fmt, & & d_cp_coo_to_coo, d_cp_coo_from_coo, & & d_cp_coo_to_fmt, d_cp_coo_from_fmt, & - & d_coo_scals, d_coo_scal, d_coo_csgetrow + & d_coo_scals, d_coo_scal, d_coo_csgetrow, d_coo_sizeof interface @@ -838,6 +843,265 @@ contains end subroutine d_base_cssv + subroutine d_cssm(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + use psb_string_mod + implicit none + class(psbn_d_base_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_dpk_), intent(in), optional :: d(:) + + real(psb_dpk_), allocatable :: tmp(:,:) + Integer :: err_act, nar,nac,nc, i + character(len=1) :: side_ + character(len=20) :: name='d_cssm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + nc = min(size(x,2), size(y,2)) + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(d)) then + if (present(side)) then + side_ = side + else + side_ = 'L' + end if + + if (psb_toupper(side_) == 'R') then + if (size(d,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nac,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nac,nc),stat=info) + if (info /= 0) info = 4000 + if (info == 0) then + do i=1, nac + tmp(i,1:nc) = d(i)*x(i,1:nc) + end do + end if + if (info == 0)& + & call a%base_cssm(alpha,tmp,beta,y,info,trans) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else if (psb_toupper(side_) == 'L') then + + if (size(d,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nar,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nar,nc),stat=info) + if (info /= 0) info = 4000 + if (info == 0)& + & call a%base_cssm(done,x,dzero,tmp,info,trans) + + if (info == 0)then + do i=1, nar + tmp(i,1:nc) = d(i)*tmp(i,1:nc) + end do + end if + if (info == 0)& + & call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else + info = 31 + call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_) + goto 9999 + end if + else + ! Side is ignored in this case + call a%base_cssm(alpha,x,beta,y,info,trans) + end if + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name, a_err='base_cssm') + goto 9999 + end if + + + return + 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_cssm + + subroutine d_cssv(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + use psb_string_mod + implicit none + class(psbn_d_base_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_dpk_), intent(in), optional :: d(:) + + real(psb_dpk_), allocatable :: tmp(:) + Integer :: err_act, nar,nac,nc, i + character(len=1) :: side_ + character(len=20) :: name='d_cssm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + nc = 1 + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(d)) then + if (present(side)) then + side_ = side + else + side_ = 'L' + end if + + if (psb_toupper(side_) == 'R') then + if (size(d,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nac,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nac),stat=info) + if (info /= 0) info = 4000 + if (info == 0) tmp(1:nac) = d(1:nac)*x(1:nac) + if (info == 0)& + & call a%base_cssm(alpha,tmp,beta,y,info,trans) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else if (psb_toupper(side_) == 'L') then + if (size(d,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nar,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nar),stat=info) + if (info /= 0) info = 4000 + if (info == 0)& + & call a%base_cssm(done,x,dzero,tmp,info,trans) + + if (info == 0) tmp(1:nar) = d(1:nar)*tmp(1:nar) + if (info == 0)& + & call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else + info = 31 + call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_) + goto 9999 + end if + else + ! Side is ignored in this case + call a%base_cssm(alpha,x,beta,y,info,trans) + end if + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name, a_err='base_cssm') + goto 9999 + end if + + + return + 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_cssv + subroutine d_scals(d,a,info) use psb_error_mod @@ -908,7 +1172,7 @@ contains ! so we throw an error. info = 700 call psb_errpush(info,name,a_err=a%get_fmt()) - write(0,*) 'Got into error path',err_act,psb_act_ret_ + if (err_act /= psb_act_ret_) then call psb_error() end if @@ -961,6 +1225,18 @@ contains !==================================== + + function d_coo_sizeof(a) result(res) + implicit none + class(psbn_d_coo_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + 1 + res = res + psb_sizeof_dp * size(a%val) + res = res + psb_sizeof_int * size(a%ia) + res = res + psb_sizeof_int * size(a%ja) + + end function d_coo_sizeof + function d_coo_get_fmt(a) result(res) implicit none @@ -1597,7 +1873,6 @@ contains call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) goto 9999 endif - if (info == 0) call psb_realloc(nz_,a%ia,info) if (info == 0) call psb_realloc(nz_,a%ja,info) if (info == 0) call psb_realloc(nz_,a%val,info) @@ -1721,7 +1996,7 @@ contains character, optional, intent(in) :: trans character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc + integer :: i,j,k,m,n, nnz, ir, jc, nac, nar real(psb_dpk_) :: acc logical :: tra Integer :: err_act @@ -1735,7 +2010,19 @@ contains call psb_errpush(info,name) goto 9999 endif - + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) @@ -1765,7 +2052,7 @@ contains character, optional, intent(in) :: trans character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc, nc + integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac real(psb_dpk_), allocatable :: acc(:) logical :: tra Integer :: err_act @@ -1775,7 +2062,24 @@ contains call psb_erractionsave(err_act) - + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) if (info /= 0) goto 9999 @@ -1805,7 +2109,7 @@ contains character, optional, intent(in) :: trans character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc + integer :: i,j,k,m,n, nnz, ir, jc, nar, nac real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:) logical :: tra @@ -1821,9 +2125,21 @@ contains goto 9999 endif + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 @@ -1859,7 +2175,7 @@ contains character, optional, intent(in) :: trans character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc, nc + integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:,:) logical :: tra @@ -1875,9 +2191,21 @@ contains goto 9999 endif + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 diff --git a/base/newserial/psbn_d_coo_impl.f03 b/base/newserial/psbn_d_coo_impl.f03 index 86c30b31..0b7156c6 100644 --- a/base/newserial/psbn_d_coo_impl.f03 +++ b/base/newserial/psbn_d_coo_impl.f03 @@ -29,7 +29,6 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 @@ -48,20 +47,20 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) if (alpha == dzero) then if (beta == dzero) then do i = 1, m - y(i,:) = dzero + y(i,1:nc) = dzero enddo else do i = 1, m - y(i,:) = beta*y(i,:) + y(i,1:nc) = beta*y(i,1:nc) end do endif return end if if (beta == dzero) then - call inner_coosm(tra,a,x,y,info) + call inner_coosm(tra,a,x(:,1:nc),y(:,1:nc),info) do i = 1, m - y(i,:) = alpha*y(i,:) + y(i,1:nc) = alpha*y(i,1:nc) end do else allocate(tmp(m,nc), stat=info) @@ -71,10 +70,10 @@ subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) goto 9999 end if - tmp(1:m,:) = x(1:m,:) - call inner_coosm(tra,a,tmp,y,info) + tmp(1:m,1:nc) = x(1:m,1:nc) + call inner_coosm(tra,a,tmp(:,1:nc),y(:,1:nc),info) do i = 1, m - y(i,:) = alpha*tmp(i,:) + beta*y(i,:) + y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) end do end if @@ -135,7 +134,7 @@ contains do if (j > nnz) exit if (a%ia(j) > i) exit - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),:) j = j + 1 end do y(i,:) = x(i,:) - acc @@ -152,7 +151,7 @@ contains j = j + 1 exit end if - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),:) j = j + 1 end do end do @@ -185,7 +184,7 @@ contains j = j - 1 exit end if - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),:) j = j - 1 end do end do @@ -396,7 +395,7 @@ contains do if (j > nnz) exit if (a%ia(j) > i) exit - acc = acc + a%val(j)*x(a%ja(j)) + acc = acc + a%val(j)*y(a%ja(j)) j = j + 1 end do y(i) = x(i) - acc @@ -413,7 +412,7 @@ contains j = j + 1 exit end if - acc = acc + a%val(j)*x(a%ja(j)) + acc = acc + a%val(j)*y(a%ja(j)) j = j + 1 end do end do @@ -427,7 +426,7 @@ contains do if (j < 1) exit if (a%ia(j) < i) exit - acc = acc + a%val(j)*x(a%ja(j)) + acc = acc + a%val(j)*y(a%ja(j)) j = j - 1 end do y(i) = x(i) - acc @@ -446,7 +445,7 @@ contains j = j - 1 exit end if - acc = acc + a%val(j)*x(a%ja(j)) + acc = acc + a%val(j)*y(a%ja(j)) j = j - 1 end do end do @@ -747,11 +746,11 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) if (alpha == dzero) then if (beta == dzero) then do i = 1, m - y(i,:) = dzero + y(i,1:nc) = dzero enddo else do i = 1, m - y(i,:) = beta*y(i,:) + y(i,1:nc) = beta*y(i,1:nc) end do endif return @@ -759,28 +758,28 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) if (.not.a%is_unit()) then if (beta == dzero) then do i = 1, m - y(i,:) = dzero + y(i,1:nc) = dzero enddo else do i = 1, m - y(i,:) = beta*y(i,:) + y(i,1:nc) = beta*y(i,1:nc) end do endif else if (a%is_unit()) then if (beta == dzero) then do i = 1, min(m,n) - y(i,:) = alpha*x(i,:) + y(i,1:nc) = alpha*x(i,1:nc) enddo do i = min(m,n)+1, m - y(i,:) = dzero + y(i,1:nc) = dzero enddo else do i = 1, min(m,n) - y(i,:) = beta*y(i,:) + alpha*x(i,:) + y(i,1:nc) = beta*y(i,1:nc) + alpha*x(i,1:nc) end do do i = min(m,n)+1, m - y(i,:) = beta*y(i,:) + y(i,1:nc) = beta*y(i,1:nc) enddo endif @@ -796,15 +795,15 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) acc = dzero do if (i>nnz) then - y(ir,:) = y(ir,:) + alpha * acc + y(ir,1:nc) = y(ir,1:nc) + alpha * acc exit endif if (a%ia(i) /= ir) then - y(ir,:) = y(ir,:) + alpha * acc + y(ir,1:nc) = y(ir,1:nc) + alpha * acc ir = a%ia(i) acc = dzero endif - acc = acc + a%val(i) * x(a%ja(i),:) + acc = acc + a%val(i) * x(a%ja(i),1:nc) i = i + 1 enddo end if @@ -815,7 +814,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) do i=1,nnz ir = a%ja(i) jc = a%ia(i) - y(ir,:) = y(ir,:) + a%val(i)*x(jc,:) + y(ir,1:nc) = y(ir,1:nc) + a%val(i)*x(jc,1:nc) enddo else if (alpha.eq.-done) then @@ -823,7 +822,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) do i=1,nnz ir = a%ja(i) jc = a%ia(i) - y(ir,:) = y(ir,:) - a%val(i)*x(jc,:) + y(ir,1:nc) = y(ir,1:nc) - a%val(i)*x(jc,1:nc) enddo else @@ -831,7 +830,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) do i=1,nnz ir = a%ja(i) jc = a%ia(i) - y(ir,:) = y(ir,:) + alpha*a%val(i)*x(jc,:) + y(ir,1:nc) = y(ir,1:nc) + alpha*a%val(i)*x(jc,1:nc) enddo end if !.....end testing on alpha @@ -861,7 +860,6 @@ function d_coo_csnmi_impl(a) result(res) integer :: i,j,k,m,n, nnz, ir, jc, nc real(psb_dpk_) :: acc - real(psb_dpk_), allocatable :: tmp(:,:) logical :: tra Integer :: err_act character(len=20) :: name='d_base_csnmi' @@ -1028,7 +1026,6 @@ contains irw = imin lrw = imax if (irw<0) then - write(debug_unit,*) ' spgtrow Error : idx no good ',irw info = 2 return end if @@ -1237,13 +1234,11 @@ subroutine d_coo_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) nza = a%get_nzeros() isza = a%get_size() -!!$ write(0,*) 'On entry to csput_impl: ',nza if (a%is_bld()) then ! Build phase. Must handle reallocations in a sensible way. if (isza < (nza+nz)) then call a%reallocate(max(nza+nz,int(1.5*isza))) isza = a%get_size() -!!$ write(0,*) 'isza: ',isza,nza+nz endif call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,& @@ -1327,7 +1322,6 @@ contains if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then nza = nza + 1 if (nza > maxsz) then - write(0,*) 'err -92 ',nza,maxsz info = -92 return endif diff --git a/base/newserial/psbn_d_csr_impl.f03 b/base/newserial/psbn_d_csr_impl.f03 index 3e2a213a..e1e0936f 100644 --- a/base/newserial/psbn_d_csr_impl.f03 +++ b/base/newserial/psbn_d_csr_impl.f03 @@ -40,7 +40,6 @@ subroutine d_csr_csmv_impl(alpha,a,x,beta,y,info,trans) end if if (.not.a%is_asb()) then - write(0,*) 'Error: csmv called on an unassembled mat' info = 1121 call psb_errpush(info,name) goto 9999 @@ -567,7 +566,6 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) info = 0 call psb_erractionsave(err_act) - if (present(trans)) then trans_ = trans else @@ -588,7 +586,7 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) goto 9999 end if - + if (alpha == dzero) then if (beta == dzero) then do i = 1, m @@ -603,18 +601,29 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) end if if (beta == dzero) then - call inner_csrsv(tra,a,x,y) - do i = 1, m - y(i) = alpha*y(i) - end do +!!$ call inner_csrsv(tra,a,x,y) + call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& + & a%irp,a%ja,a%val,x,y) + if (alpha == done) then + ! do nothing + else if (alpha == -done) then + do i = 1, m + y(i) = -y(i) + end do + else + do i = 1, m + y(i) = alpha*y(i) + end do + end if else allocate(tmp(m), stat=info) if (info /= 0) then - write(0,*) 'Memory allocation error in CSRSV ' return end if tmp(1:m) = x(1:m) - call inner_csrsv(tra,a,tmp,y) + call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& + & a%irp,a%ja,a%val,tmp,y) +!!$ call inner_csrsv(tra,a,tmp,y) do i = 1, m y(i) = alpha*tmp(i) + beta*y(i) end do @@ -634,53 +643,56 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans) contains - subroutine inner_csrsv(tra,a,x,y) + subroutine inner_csrsv(tra,lower,unit,n,irp,ja,val,x,y) implicit none - logical, intent(in) :: tra - class(psbn_d_csr_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:) - real(psb_dpk_), intent(out) :: y(:) + logical, intent(in) :: tra,lower,unit +!!$ class(psbn_d_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: irp(*), ja(*),n + + real(psb_dpk_), intent(in) :: val(*) + real(psb_dpk_), intent(in) :: x(*) + real(psb_dpk_), intent(out) :: y(*) integer :: i,j,k,m, ir, jc real(psb_dpk_) :: acc if (.not.tra) then - if (a%is_lower()) then - if (a%is_unit()) then - do i=1, a%get_nrows() + if (lower) then + if (unit) then + do i=1, n acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) end do y(i) = x(i) - acc end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() + else if (.not.unit) then + do i=1, n acc = dzero - do j=a%irp(i), a%irp(i+1)-2 - acc = acc + a%val(j)*x(a%ja(j)) + do j=irp(i), irp(i+1)-2 + acc = acc + val(j)*y(ja(j)) end do - y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1) + y(i) = (x(i) - acc)/val(irp(i+1)-1) end do end if - else if (a%is_upper()) then + else if (.not.lower) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 + if (unit) then + do i=n, 1, -1 acc = dzero - do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) end do y(i) = x(i) - acc end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 + else if (.not.unit) then + do i=n, 1, -1 acc = dzero - do j=a%irp(i)+1, a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j)) + do j=irp(i)+1, irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) end do - y(i) = (x(i) - acc)/a%val(a%irp(i)) + y(i) = (x(i) - acc)/val(irp(i)) end do end if @@ -688,46 +700,46 @@ contains else if (tra) then - do i=1, a%get_nrows() + do i=1, n y(i) = x(i) end do - if (a%is_lower()) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 + if (lower) then + if (unit) then + do i=n, 1, -1 acc = y(i) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc end do end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - y(i) = y(i)/a%val(a%irp(i+1)-1) + else if (.not.unit) then + do i=n, 1, -1 + y(i) = y(i)/val(irp(i+1)-1) acc = y(i) - do j=a%irp(i), a%irp(i+1)-2 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc + do j=irp(i), irp(i+1)-2 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc end do end do end if - else if (a%is_upper()) then + else if (.not.lower) then - if (a%is_unit()) then - do i=1, a%get_nrows() + if (unit) then + do i=1, n acc = y(i) - do j=a%irp(i), a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc end do end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - y(i) = y(i)/a%val(a%irp(i)) + else if (.not.unit) then + do i=1, n + y(i) = y(i)/val(irp(i)) acc = y(i) - do j=a%irp(i)+1, a%irp(i+1)-1 - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc + do j=irp(i)+1, irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc end do end do end if @@ -779,7 +791,6 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) nc = min(size(x,2) , size(y,2)) if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 @@ -800,9 +811,11 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) end if if (beta == dzero) then - call inner_csrsm(tra,a,x,y,info) +!!$ call inner_csrsm(tra,a,x,y,info) + call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1),y,size(y,1),info) do i = 1, m - y(i,:) = alpha*y(i,:) + y(i,1:nc) = alpha*y(i,1:nc) end do else allocate(tmp(m,nc), stat=info) @@ -812,10 +825,12 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) goto 9999 end if - tmp(1:m,:) = x(1:m,:) - call inner_csrsm(tra,a,tmp,y,info) + tmp(1:m,:) = x(1:m,1:nc) + call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,tmp,size(tmp,1),y,size(y,1),info) +!!$ call inner_csrsm(tra,a,tmp,y,info) do i = 1, m - y(i,:) = alpha*tmp(i,:) + beta*y(i,:) + y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) end do end if @@ -841,18 +856,19 @@ subroutine d_csr_cssm_impl(alpha,a,x,beta,y,info,trans) contains - subroutine inner_csrsm(tra,a,x,y,info) + subroutine inner_csrsm(tra,lower,unit,nr,nc,& + & irp,ja,val,x,ldx,y,ldy,info) implicit none - logical, intent(in) :: tra - class(psbn_d_csr_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:,:) - real(psb_dpk_), intent(out) :: y(:,:) + logical, intent(in) :: tra,lower,unit + integer, intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*) + real(psb_dpk_), intent(in) :: val(*), x(ldx,*) + real(psb_dpk_), intent(out) :: y(ldy,*) integer, intent(out) :: info integer :: i,j,k,m, ir, jc real(psb_dpk_), allocatable :: acc(:) info = 0 - allocate(acc(size(x,2)), stat=info) + allocate(acc(nc), stat=info) if(info /= 0) then info=4010 return @@ -861,41 +877,41 @@ contains if (.not.tra) then - if (a%is_lower()) then - if (a%is_unit()) then - do i=1, a%get_nrows() + if (lower) then + if (unit) then + do i=1, nr acc = dzero do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),1:nc) end do - y(i,:) = x(i,:) - acc + y(i,1:nc) = x(i,1:nc) - acc end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() + else if (.not.unit) then + do i=1, nr acc = dzero do j=a%irp(i), a%irp(i+1)-2 - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),1:nc) end do - y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1) + y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i+1)-1) end do end if - else if (a%is_upper()) then + else if (.not.lower) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 + if (unit) then + do i=nr, 1, -1 acc = dzero do j=a%irp(i), a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),1:nc) end do - y(i,:) = x(i,:) - acc + y(i,1:nc) = x(i,1:nc) - acc end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 + else if (.not.unit) then + do i=nr, 1, -1 acc = dzero do j=a%irp(i)+1, a%irp(i+1)-1 - acc = acc + a%val(j)*x(a%ja(j),:) + acc = acc + a%val(j)*y(a%ja(j),1:nc) end do - y(i,:) = (x(i,:) - acc)/a%val(a%irp(i)) + y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i)) end do end if @@ -903,46 +919,46 @@ contains else if (tra) then - do i=1, a%get_nrows() - y(i,:) = x(i,:) + do i=1, nr + y(i,1:nc) = x(i,1:nc) end do - if (a%is_lower()) then - if (a%is_unit()) then - do i=a%get_nrows(), 1, -1 - acc = y(i,:) + if (lower) then + if (unit) then + do i=nr, 1, -1 + acc = y(i,1:nc) do j=a%irp(i), a%irp(i+1)-1 jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc end do end do - else if (.not.a%is_unit()) then - do i=a%get_nrows(), 1, -1 - y(i,:) = y(i,:)/a%val(a%irp(i+1)-1) - acc = y(i,:) + else if (.not.unit) then + do i=nr, 1, -1 + y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i+1)-1) + acc = y(i,1:nc) do j=a%irp(i), a%irp(i+1)-2 jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc end do end do end if - else if (a%is_upper()) then + else if (.not.lower) then - if (a%is_unit()) then - do i=1, a%get_nrows() - acc = y(i,:) + if (unit) then + do i=1, nr + acc = y(i,1:nc) do j=a%irp(i), a%irp(i+1)-1 jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc end do end do - else if (.not.a%is_unit()) then - do i=1, a%get_nrows() - y(i,:) = y(i,:)/a%val(a%irp(i)) - acc = y(i,:) + else if (.not.unit) then + do i=1, nr + y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i)) + acc = y(i,1:nc) do j=a%irp(i)+1, a%irp(i+1)-1 jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc end do end do end if @@ -1120,7 +1136,6 @@ contains irw = imin lrw = min(imax,a%get_nrows()) if (irw<0) then - write(debug_unit,*) ' spgtrow Error : idx no good ',irw info = 2 return end if @@ -1131,13 +1146,14 @@ contains nzin_ = 0 endif - nzt = a%irp(lrw)-a%irp(irw) + nzt = a%irp(lrw+1)-a%irp(irw) nz = 0 - + call psb_ensure_size(nzin_+nzt,ia,info) if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return if (present(iren)) then @@ -1671,6 +1687,9 @@ subroutine d_mv_csr_to_fmt_impl(a,b,info) select type (b) class is (psbn_d_coo_sparse_mat) call a%mv_to_coo(b,info) + ! Need to fix trivial copies! +!!$ class is (psbn_d_csr_sparse_mat) +!!$ call a%mv_to_coo(b,info) class default call tmp%mv_from_fmt(a,info) if (info == 0) call b%mv_from_coo(tmp,info) diff --git a/base/newserial/psbn_d_csr_mat_mod.f03 b/base/newserial/psbn_d_csr_mat_mod.f03 index c3a179e5..50e00e47 100644 --- a/base/newserial/psbn_d_csr_mat_mod.f03 +++ b/base/newserial/psbn_d_csr_mat_mod.f03 @@ -34,6 +34,7 @@ module psbn_d_csr_mat_mod 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 end type psbn_d_csr_sparse_mat private :: d_csr_get_nzeros, d_csr_csmm, d_csr_csmv, d_csr_cssm, d_csr_cssv, & & d_csr_csput, d_csr_reallocate_nz, d_csr_allocate_mnnz, & @@ -42,7 +43,8 @@ module psbn_d_csr_mat_mod & d_mv_csr_to_coo, d_mv_csr_from_coo, & & d_cp_csr_to_fmt, d_cp_csr_from_fmt, & & d_mv_csr_to_fmt, d_mv_csr_from_fmt, & - & d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_get_size + & d_csr_scals, d_csr_scal, d_csr_trim, d_csr_csgetrow, d_csr_get_size, & + & d_csr_sizeof interface @@ -234,6 +236,18 @@ contains ! !===================================== + + function d_csr_sizeof(a) result(res) + implicit none + class(psbn_d_csr_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + psb_sizeof_dp * size(a%val) + res = res + psb_sizeof_int * size(a%irp) + res = res + psb_sizeof_int * size(a%ja) + + end function d_csr_sizeof + function d_csr_get_fmt(a) result(res) implicit none class(psbn_d_csr_sparse_mat), intent(in) :: a @@ -1151,7 +1165,6 @@ contains if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 @@ -1205,7 +1218,6 @@ contains if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' info = 1121 call psb_errpush(info,name) goto 9999 diff --git a/base/newserial/psbn_mat_mod.f03 b/base/newserial/psbn_mat_mod.f03 index 3c5b18a1..a6da6f32 100644 --- a/base/newserial/psbn_mat_mod.f03 +++ b/base/newserial/psbn_mat_mod.f03 @@ -40,6 +40,8 @@ module psbn_d_mat_mod procedure, pass(a) :: is_triangle procedure, pass(a) :: is_unit procedure, pass(a) :: get_fmt => sparse_get_fmt + procedure, pass(a) :: sizeof => d_sizeof + ! Memory/data management procedure, pass(a) :: csall @@ -78,12 +80,32 @@ module psbn_d_mat_mod & is_unit, get_neigh, csall, csput, d_csgetrow,& & d_csgetblk, csclip, d_cscnv, d_cscnv_ip, & & reallocate_nz, free, trim, & - & d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, & + & sparse_print, & & set_nrows, set_ncols, set_dupl, & & set_state, set_null, set_bld, & & set_upd, set_asb, set_sorted, & & set_upper, set_lower, set_triangle, & - & set_unit, csnmi, get_diag, d_scals, d_scal + & set_unit, get_diag + + interface psb_sizeof + module procedure d_sizeof + end interface + + interface psbn_csmm + module procedure d_csmm, d_csmv + end interface + + interface psbn_cssm + module procedure d_cssm, d_cssv + end interface + + interface psbn_csnmi + module procedure csnmi + end interface + + interface psbn_scal + module procedure d_scals, d_scal + end interface contains @@ -100,6 +122,20 @@ contains ! !===================================== + + function d_sizeof(a) result(res) + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + + res = 0 + if (allocated(a%a)) then + res = a%a%sizeof() + end if + + end function d_sizeof + + function sparse_get_fmt(a) result(res) implicit none @@ -1265,7 +1301,7 @@ contains call move_alloc(altmp,b%a) call b%set_asb() - + call b%trim() call psb_erractionrestore(err_act) return @@ -1357,7 +1393,7 @@ contains call move_alloc(altmp,a%a) call a%set_asb() - + call a%trim() call psb_erractionrestore(err_act) return @@ -1460,14 +1496,15 @@ contains end subroutine d_csmv - subroutine d_cssm(alpha,a,x,beta,y,info,trans) + subroutine d_cssm(alpha,a,x,beta,y,info,trans,side,d) use psb_error_mod implicit none class(psbn_d_sparse_mat), intent(in) :: a real(kind(1.d0)), intent(in) :: alpha, beta, x(:,:) real(kind(1.d0)), intent(inout) :: y(:,:) integer, intent(out) :: info - character, optional, intent(in) :: trans + character, optional, intent(in) :: trans, side + real(psb_dpk_), intent(in), optional :: d(:) Integer :: err_act character(len=20) :: name='psbn_cssm' logical, parameter :: debug=.false. @@ -1480,7 +1517,7 @@ contains goto 9999 endif - call a%a%cssm(alpha,x,beta,y,info,trans) + call a%a%cssm(alpha,x,beta,y,info,trans,side,d) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -1497,14 +1534,15 @@ contains end subroutine d_cssm - subroutine d_cssv(alpha,a,x,beta,y,info,trans) + subroutine d_cssv(alpha,a,x,beta,y,info,trans,side,d) use psb_error_mod implicit none class(psbn_d_sparse_mat), intent(in) :: a real(kind(1.d0)), intent(in) :: alpha, beta, x(:) real(kind(1.d0)), intent(inout) :: y(:) integer, intent(out) :: info - character, optional, intent(in) :: trans + character, optional, intent(in) :: trans, side + real(psb_dpk_), intent(in), optional :: d(:) Integer :: err_act character(len=20) :: name='psbn_cssv' logical, parameter :: debug=.false. @@ -1516,8 +1554,8 @@ contains call psb_errpush(info,name) goto 9999 endif - - call a%a%cssm(alpha,x,beta,y,info,trans) + + call a%a%cssm(alpha,x,beta,y,info,trans,side,d) if (info /= 0) goto 9999 diff --git a/base/psblas/psb_dnrmi.f90 b/base/psblas/psb_dnrmi.f90 index 32b40691..7cbd515e 100644 --- a/base/psblas/psb_dnrmi.f90 +++ b/base/psblas/psb_dnrmi.f90 @@ -95,7 +95,7 @@ function psb_dnrmi(a,desc_a,info) end if if ((m /= 0).and.(n /= 0)) then - nrmi = a%csnmi() + nrmi = psbn_csnmi(a) if(info /= 0) then info=4010 ch_err='psb_csnmi' diff --git a/base/psblas/psb_dspmm.f90 b/base/psblas/psb_dspmm.f90 index 9737e501..cda9d00d 100644 --- a/base/psblas/psb_dspmm.f90 +++ b/base/psblas/psb_dspmm.f90 @@ -251,7 +251,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& if(info /= 0) exit blk ! local Matrix-vector product - call a%csmm(alpha,x(:,jjx+i-1:jjx+i-1+ib-1),& + call psbn_csmm(alpha,a,x(:,jjx+i-1:jjx+i-1+ib-1),& & beta,y(:,jjy+i-1:jjy+i-1+ib-1),info,trans=trans_) if(info /= 0) exit blk @@ -266,7 +266,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& if (doswap_)& & call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& & ib1,dzero,x(:,1:ik),desc_a,iwork,info) - if (info == 0) call a%csmm(alpha,x(:,1:ik),beta,y(:,1:ik),info) + if (info == 0) call psbn_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info) end if if(info /= 0) then info = 4011 @@ -313,7 +313,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,& if (info == 0) call psi_ovrl_upd(x,desc_a,psb_avg_,info) y(nrow+1:ncol,1:ik) = dzero - if (info == 0) call a%csmm(alpha,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) + if (info == 0) call psbn_csmm(alpha,a,x(:,1:ik),beta,y(:,1:ik),info,trans=trans_) if (debug_level >= psb_debug_comp_) & & write(debug_unit,*) me,' ',trim(name),' csmm ', info if (info /= 0) then @@ -584,7 +584,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& & dzero,x,desc_a,iwork,info,data=psb_comm_halo_) end if ! Just for fun - call a%csmm(alpha,x,beta,y,info) + call psbn_csmm(alpha,a,x,beta,y,info) if(info /= 0) then info = 4011 @@ -633,7 +633,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,& yp(nrow+1:ncol) = dzero ! local Matrix-vector product - if (info == 0) call a%csmm(alpha,x,beta,y,info,trans=trans_) + if (info == 0) call psbn_csmm(alpha,a,x,beta,y,info,trans=trans_) if (debug_level >= psb_debug_comp_) & & write(debug_unit,*) me,' ',trim(name),' csmm ', info diff --git a/base/psblas/psb_dspsm.f90 b/base/psblas/psb_dspsm.f90 index 6c7cb8d2..fafe4773 100644 --- a/base/psblas/psb_dspsm.f90 +++ b/base/psblas/psb_dspsm.f90 @@ -64,7 +64,7 @@ ! desc_a - type(psb_desc_type). The communication descriptor. ! info - integer. Return code ! trans - character(optional). Whether A or A'. If not present 'N' is assumed. -! unitd - character(optional). Specify some type of operation with +! side - character(optional). Specify some type of operation with ! the diagonal matrix D. ! choice - integer(optional). The kind of update to perform on overlap elements. ! d(:) - real , optional Matrix for diagonal scaling. @@ -75,7 +75,7 @@ ! ! subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& - & trans, unitd, choice, diag, k, jx, jy, work) + & trans, side, choice, diag, k, jx, jy, work) use psb_spmat_type use psb_serial_mod @@ -86,17 +86,18 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& use psb_error_mod use psb_string_mod use psb_penv_mod + use psbn_d_mat_mod implicit none real(psb_dpk_), intent(in) :: alpha, beta real(psb_dpk_), intent(in), target :: x(:,:) real(psb_dpk_), intent(inout), target :: y(:,:) - type (psb_dspmat_type), intent(in) :: a + type(psbn_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(psb_dpk_), intent(in), optional, target :: diag(:) real(psb_dpk_), optional, target :: work(:) - character, intent(in), optional :: trans, unitd + character, intent(in), optional :: trans, side integer, intent(in), optional :: choice integer, intent(in), optional :: k, jx, jy @@ -106,7 +107,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& & ix, iy, ik, ijx, ijy, i, lld,& & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm - character :: lunitd + character :: lside integer, parameter :: nb=4 real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:) character :: itrans @@ -158,10 +159,10 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& choice_ = psb_avg_ endif - if (present(unitd)) then - lunitd = psb_toupper(unitd) + if (present(side)) then + lside = psb_toupper(side) else - lunitd = 'U' + lside = 'U' endif if (present(trans)) then @@ -192,8 +193,6 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& ! check for presence/size of a work area iwork => null() liwork= 2*ncol - if (a%pr(1) /= 0) llwork = liwork + m * ik - if (a%pl(1) /= 0) llwork = llwork + m * ik if (present(work)) then if (size(work) >= liwork) then aliw =.false. @@ -259,7 +258,7 @@ subroutine psb_dspsm(alpha,a,x,beta,y,desc_a,info,& ! Perform local triangular system solve xp => x(iix:lldx,jjx:jjx+ik-1) yp => y(iiy:lldy,jjy:jjy+ik-1) - call a%cssm(alpha,xp,beta,yp,info,unitd=lunitd,d=id,trans=itrans) + call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=diag,trans=trans) if(info /= 0) then info = 4010 @@ -357,14 +356,14 @@ end subroutine psb_dspsm ! desc_a - type(psb_desc_type). The communication descriptor. ! info - integer. Return code ! trans - character(optional). Whether A or A'. If not present 'N' is assumed. -! unitd - character(optional). Specify some type of operation with +! side - character(optional). Specify some type of operation with ! the diagonal matrix D. ! choice - integer(optional). The kind of update to perform on overlap elements. ! d(:) - real , optional Matrix for diagonal scaling. ! work(:) - real , optional Working area. ! subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& - & trans, unitd, choice, diag, work) + & trans, side, choice, diag, work) use psb_spmat_type use psb_serial_mod use psb_descriptor_type @@ -374,17 +373,18 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& use psb_error_mod use psb_string_mod use psb_penv_mod + use psbn_d_mat_mod implicit none real(psb_dpk_), intent(in) :: alpha, beta real(psb_dpk_), intent(in), target :: x(:) real(psb_dpk_), intent(inout), target :: y(:) - type(psb_dspmat_type), intent(in) :: a + type(psbn_d_sparse_mat), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info real(psb_dpk_), intent(in), optional, target :: diag(:) real(psb_dpk_), optional, target :: work(:) - character, intent(in), optional :: trans, unitd + character, intent(in), optional :: trans, side integer, intent(in), optional :: choice ! locals @@ -393,7 +393,7 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& & ix, iy, ik, jx, jy, i, lld,& & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm - character :: lunitd + character :: lside integer, parameter :: nb=4 real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:), id(:) character :: itrans @@ -429,10 +429,10 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& choice_ = psb_avg_ endif - if (present(unitd)) then - lunitd = psb_toupper(unitd) + if (present(side)) then + lside = psb_toupper(side) else - lunitd = 'U' + lside = 'U' endif if (present(trans)) then @@ -529,7 +529,8 @@ subroutine psb_dspsv(alpha,a,x,beta,y,desc_a,info,& ! Perform local triangular system solve xp => x(iix:lldx) yp => y(iiy:lldy) - call a%cssm(alpha,xp,beta,yp,info,unitd=lunitd,d=id,trans=itrans) + call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=diag,trans=trans) +!!$ call psbn_cssm(alpha,a,xp,beta,yp,info,side=side,d=id,trans=itrans) if(info /= 0) then info = 4010 diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 5f128885..a89af367 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -109,7 +109,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) & write(debug_unit,*) me,' ',trim(name),':allocating size:',length_ia1 !....allocate aspk, ia1, ia2..... - call a%csall(loc_row,loc_col,length_ia1,info) + call a%csall(loc_row,loc_col,info,nz=length_ia1) if(info /= 0) then info=4010 ch_err='sp_all' diff --git a/prec/psb_dbjac_aply.f90 b/prec/psb_dbjac_aply.f90 index 23877eed..c4502946 100644 --- a/prec/psb_dbjac_aply.f90 +++ b/prec/psb_dbjac_aply.f90 @@ -38,6 +38,7 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod + use psbn_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 a9973691..5566fdbf 100644 --- a/prec/psb_dbjac_bld.f90 +++ b/prec/psb_dbjac_bld.f90 @@ -30,6 +30,7 @@ !!$ !!$ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) + use psbn_d_mat_mod use psb_base_mod use psb_prec_mod, psb_protect_name => psb_dbjac_bld implicit none @@ -37,7 +38,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) ! .. Scalar Arguments .. integer, intent(out) :: info ! .. array Arguments .. - type(psb_dspmat_type), intent(in), target :: a + type(psbn_d_sparse_mat), intent(in), target :: a type(psb_dprec_type), intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd @@ -47,6 +48,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) integer :: int_err(5) character :: trans, unitd type(psb_dspmat_type) :: atmp + type(psbn_d_csr_sparse_mat), allocatable :: lf, uf real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8 integer nztota, err_act, n_row, nrow_a,n_col, nhalo integer :: ictxt,np,me @@ -61,7 +63,7 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) - m = a%m + m = a%get_nrows() if (m < 0) then info = 10 int_err(1) = 1 @@ -89,12 +91,12 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) if (allocated(p%av)) then if (size(p%av) < psb_bp_ilu_avsz) then do i=1,size(p%av) - call psb_sp_free(p%av(i),info) - if (info /= 0) then - ! Actually, we don't care here about this. - ! Just let it go. - ! return - end if + call p%av(i)%free() +!!$ if (info /= 0) then +!!$ ! Actually, we don't care here about this. +!!$ ! Just let it go. +!!$ ! return +!!$ end if enddo deallocate(p%av,stat=info) endif @@ -108,17 +110,19 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) endif nrow_a = psb_cd_get_local_rows(desc_a) - nztota = psb_sp_get_nnzeros(a) + nztota = a%get_nzeros() n_col = psb_cd_get_local_cols(desc_a) nhalo = n_col-nrow_a n_row = p%desc_data%matrix_data(psb_n_row_) - p%av(psb_l_pr_)%m = n_row - p%av(psb_l_pr_)%k = n_row - p%av(psb_u_pr_)%m = n_row - p%av(psb_u_pr_)%k = n_row - call psb_sp_all(n_row,n_row,p%av(psb_l_pr_),nztota,info) - if (info == 0) call psb_sp_all(n_row,n_row,p%av(psb_u_pr_),nztota,info) + + allocate(lf,uf,stat=info) + if (info == 0) call lf%allocate(n_row,n_row,nztota) + if (info == 0) call uf%allocate(n_row,n_row,nztota) + + +!!$ call p%av(psb_l_pr_)%csall(n_row,n_row,info,nztota) +!!$ if (info == 0) call p%av(psb_u_pr_)%csall(n_row,n_row,info,nztota) if(info/=0) then info=4010 ch_err='psb_sp_all' @@ -140,24 +144,34 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info) endif t3 = psb_wtime() - ! This is where we have mo renumbering, thus no need + ! This is where we have no renumbering, thus no need ! for ATMP - - call psb_ilu_fct(a,p%av(psb_l_pr_),p%av(psb_u_pr_),p%d,info) - if(info/=0) then +!!$ call p%av(psb_l_pr_)%cscnv(info,type='CSR') +!!$ call p%av(psb_u_pr_)%cscnv(info,type='CSR') + + call psb_ilu_fct(a,lf,uf,p%d,info) + + if(info==0) then + call move_alloc(lf,p%av(psb_l_pr_)%a) + call move_alloc(uf,p%av(psb_u_pr_)%a) + call p%av(psb_l_pr_)%set_asb() + call p%av(psb_u_pr_)%set_asb() + call p%av(psb_l_pr_)%trim() + call p%av(psb_u_pr_)%trim() + else info=4010 ch_err='psb_ilu_fct' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - - if (psb_sp_getifld(psb_upd_,p%av(psb_u_pr_),info) /= psb_upd_perm_) then - call psb_sp_trim(p%av(psb_u_pr_),info) - endif - - if (psb_sp_getifld(psb_upd_,p%av(psb_l_pr_),info) /= psb_upd_perm_) then - call psb_sp_trim(p%av(psb_l_pr_),info) - endif +!!$ +!!$ if (psb_sp_getifld(psb_upd_,p%av(psb_u_pr_),info) /= psb_upd_perm_) then +!!$ call psb_sp_trim(p%av(psb_u_pr_),info) +!!$ endif +!!$ +!!$ if (psb_sp_getifld(psb_upd_,p%av(psb_l_pr_),info) /= psb_upd_perm_) then +!!$ call psb_sp_trim(p%av(psb_l_pr_),info) +!!$ endif case(psb_f_none_) diff --git a/prec/psb_dilu_fct.f90 b/prec/psb_dilu_fct.f90 index 2fd528a7..88d5043e 100644 --- a/prec/psb_dilu_fct.f90 +++ b/prec/psb_dilu_fct.f90 @@ -37,12 +37,13 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck) ! ! use psb_base_mod + use psbn_d_mat_mod implicit none ! .. Scalar Arguments .. integer, intent(out) :: info ! .. Array Arguments .. - type(psb_dspmat_type),intent(in) :: a - type(psb_dspmat_type),intent(inout) :: l,u + type(psbn_d_sparse_mat),intent(in) :: a + type(psbn_d_csr_sparse_mat),intent(inout) :: l,u type(psb_dspmat_type),intent(in), optional, target :: blck real(psb_dpk_), intent(inout) :: d(:) ! .. Local Scalars .. @@ -76,25 +77,26 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck) blck_%m=0 endif - call psb_dilu_fctint(m,a%m,a,blck_%m,blck_,& - & d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info) + call psb_dilu_fctint(m,a%get_nrows(),a,blck_%m,blck_,& + & d,l%val,l%ja,l%irp,u%val,u%ja,u%irp,l1,l2,info) if(info /= 0) then info=4010 ch_err='psb_dilu_fctint' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - - l%infoa(1) = l1 - l%fida = 'CSR' - l%descra = 'TLU' - u%infoa(1) = l2 - u%fida = 'CSR' - u%descra = 'TUU' - l%m = m - l%k = m - u%m = m - u%k = m + + call l%set_triangle() + call l%set_lower() + call l%set_unit() + call u%set_triangle() + call u%set_upper() + call u%set_unit() + call l%set_nrows(m) + call l%set_ncols(m) + call u%set_nrows(m) + call u%set_ncols(m) + if (present(blck)) then blck_ => null() else @@ -124,17 +126,22 @@ contains & d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info) implicit none - type(psb_dspmat_type) :: a,b + type(psbn_d_sparse_mat) :: a + type(psb_dspmat_type) :: b integer :: m,ma,mb,l1,l2,info integer, dimension(:) :: lia1,lia2,uia1,uia2 real(psb_dpk_), dimension(:) :: laspk,uaspk,d - integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act + integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act, nz real(psb_dpk_) :: dia,temp integer, parameter :: nrb=16 type(psb_dspmat_type) :: trw + integer, allocatable :: irow(:), icol(:) + real(psb_dpk_), allocatable :: val(:) + integer :: int_err(5) character(len=20) :: name, ch_err + name='psb_dilu_fctint' if(psb_get_errstatus() /= 0) return @@ -162,61 +169,23 @@ contains d(i) = dzero ! - ! Here we take a fast shortcut if possible, otherwise - ! use spgtblk, slower but able (in principle) to handle - ! anything. ! - if (a%fida=='CSR') then - do j = a%ia2(i), a%ia2(i+1) - 1 - k = a%ia1(j) - ! write(0,*)'KKKKK',k - if ((k < i).and.(k >= 1)) then - l1 = l1 + 1 - laspk(l1) = a%aspk(j) - lia1(l1) = k - else if (k == i) then - d(i) = a%aspk(j) - else if ((k > i).and.(k <= m)) then - l2 = l2 + 1 - uaspk(l2) = a%aspk(j) - uia1(l2) = k - end if - enddo - - else - - if ((mod(i,nrb) == 1).or.(nrb==1)) then - irb = min(ma-i+1,nrb) - call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1) - if(info /= 0) then - info=4010 - ch_err='psb_sp_getblk' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ktrw=1 + call a%csget(i,i,nz,irow,icol,val,info) + do j=1, nz + k = icol(j) + ! write(0,*)'KKKKK',k + if ((k < i).and.(k >= 1)) then + l1 = l1 + 1 + laspk(l1) = val(j) + lia1(l1) = k + else if (k == i) then + d(i) = val(j) + else if ((k > i).and.(k <= m)) then + l2 = l2 + 1 + uaspk(l2) = val(j) + uia1(l2) = k end if - - do - if (ktrw > trw%infoa(psb_nnz_)) exit - if (trw%ia1(ktrw) > i) exit - k = trw%ia2(ktrw) - if ((k < i).and.(k >= 1)) then - l1 = l1 + 1 - laspk(l1) = trw%aspk(ktrw) - lia1(l1) = k - else if (k == i) then - d(i) = trw%aspk(ktrw) - else if ((k > i).and.(k <= m)) then - l2 = l2 + 1 - uaspk(l2) = trw%aspk(ktrw) - uia1(l2) = k - end if - ktrw = ktrw + 1 - enddo - - end if - + end do !!$ lia2(i+1) = l1 + 1 diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index 68c5b711..4ab1bb18 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -115,7 +115,7 @@ subroutine psb_dprecbld(aa,desc_a,p,info,upd) call psb_check_def(p%iprcparm(psb_f_type_),'fact',& & psb_f_ilu_n_,is_legal_ml_fact) - call psb_bjac_bld(a,desc_a,p,upd_,info) + call psb_bjac_bld(aa,desc_a,p,upd_,info) if(info /= 0) then call psb_errpush(4010,name,a_err='psb_bjac_bld') diff --git a/prec/psb_dprecinit.f90 b/prec/psb_dprecinit.f90 index 513013aa..2ccfc58e 100644 --- a/prec/psb_dprecinit.f90 +++ b/prec/psb_dprecinit.f90 @@ -55,12 +55,12 @@ subroutine psb_dprecinit(p,ptype,info) p%iprcparm(psb_p_type_) = psb_diag_ p%iprcparm(psb_f_type_) = psb_f_none_ -!!$ case ('BJAC') -!!$ p%iprcparm(:) = 0 -!!$ p%iprcparm(psb_p_type_) = psb_bjac_ -!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_ -!!$ p%iprcparm(psb_ilu_fill_in_) = 0 -!!$ + case ('BJAC') + p%iprcparm(:) = 0 + p%iprcparm(psb_p_type_) = psb_bjac_ + p%iprcparm(psb_f_type_) = psb_f_ilu_n_ + p%iprcparm(psb_ilu_fill_in_) = 0 + case default write(0,*) 'Unknown preconditioner type request "',ptype,'"' info = 2 diff --git a/prec/psb_prec_mod.f90 b/prec/psb_prec_mod.f90 index d4e22b4d..0388ad76 100644 --- a/prec/psb_prec_mod.f90 +++ b/prec/psb_prec_mod.f90 @@ -330,9 +330,12 @@ module psb_prec_mod 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 psbn_d_mat_mod integer, intent(out) :: info - type(psb_dspmat_type),intent(in) :: a - type(psb_dspmat_type),intent(inout) :: l,u + type(psbn_d_sparse_mat),intent(in) :: a + type(psbn_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 real(psb_dpk_), intent(inout) :: d(:) end subroutine psb_dilu_fct @@ -365,10 +368,11 @@ module psb_prec_mod character, intent(in) :: upd end subroutine psb_sbjac_bld subroutine psb_dbjac_bld(a,desc_a,p,upd,info) + use psbn_d_mat_mod use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ use psb_prec_type, only : psb_dprec_type integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a + type(psbn_d_sparse_mat), intent(in), target :: a type(psb_dprec_type), intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd diff --git a/prec/psb_prec_type.f90 b/prec/psb_prec_type.f90 index ab4942c3..2aadfea9 100644 --- a/prec/psb_prec_type.f90 +++ b/prec/psb_prec_type.f90 @@ -41,6 +41,7 @@ module psb_prec_type & psb_dspmat_type, psb_zspmat_type, psb_dpk_, psb_spk_, psb_long_int_k_,& & psb_desc_type, psb_sizeof, psb_sp_free, psb_cdfree,& & psb_erractionsave, psb_erractionrestore, psb_error, psb_get_errstatus + use psbn_d_mat_mod, only : psbn_d_sparse_mat integer, parameter :: psb_min_prec_=0, psb_noprec_=0, psb_diag_=1, & & psb_bjac_=2, psb_max_prec_=2 @@ -75,7 +76,7 @@ module psb_prec_type end type psb_sprec_type type psb_dprec_type - type(psb_dspmat_type), allocatable :: av(:) + type(psbn_d_sparse_mat), allocatable :: av(:) real(psb_dpk_), allocatable :: d(:) type(psb_desc_type) :: desc_data integer, allocatable :: iprcparm(:) @@ -402,7 +403,8 @@ contains if (allocated(p%av)) then do i=1,size(p%av) - call psb_sp_free(p%av(i),info) +!!$ call psb_sp_free(p%av(i),info) + call p%av(i)%free() if (info /= 0) then ! Actually, we don't care here about this. ! Just let it go. @@ -600,6 +602,7 @@ contains function psb_dprec_sizeof(prec) result(val) + use psbn_d_mat_mod type(psb_dprec_type), intent(in) :: prec integer(psb_long_int_k_) :: val integer :: i diff --git a/test/pargen/ppde.f90 b/test/pargen/ppde.f90 index 29b28a53..663d9df6 100644 --- a/test/pargen/ppde.f90 +++ b/test/pargen/ppde.f90 @@ -94,7 +94,7 @@ program ppde real(psb_dpk_) :: err, eps ! other variables - integer :: info + integer :: info, i character(len=20) :: name,ch_err info=0 @@ -131,7 +131,6 @@ program ppde call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) write(*,'(" ")') ! @@ -156,7 +155,6 @@ program ppde if (iam == psb_root_) write(*,'("Preconditioner time : ",es12.5)')tprec if (iam == psb_root_) write(*,'(" ")') - ! ! iterative method parameters ! @@ -178,8 +176,7 @@ program ppde t2 = psb_wtime() - t1 call psb_amx(ictxt,t2) -!!$ amatsize = psb_sizeof(a) - amatsize = 0 + amatsize = psb_sizeof(a) descsize = psb_sizeof(desc_a) precsize = psb_sizeof(prec) call psb_sum(ictxt,amatsize) diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index 0848bcf1..6201d28c 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -1,11 +1,11 @@ 7 Number of entries below this BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES -DIAG Preconditioner NONE DIAG BJAC -CSR Storage format for matrix A: CSR COO JAD +BJAC Preconditioner NONE DIAG BJAC +COO Storage format for matrix A: CSR COO JAD 060 Domain size (acutal system is this**3) -1 Stopping criterion +2 Stopping criterion 0400 MAXIT --40 ITRACE +-01 ITRACE 20 IRST restart for RGMRES and BiCGSTABL