diff --git a/base/modules/psb_error_mod.F90 b/base/modules/psb_error_mod.F90 index 104083b7..bad1b24e 100644 --- a/base/modules/psb_error_mod.F90 +++ b/base/modules/psb_error_mod.F90 @@ -469,6 +469,8 @@ contains write (0,'("Exactly one of the optional arguments ",a," must be present")')a_e_d case(582) write (0,'("Argument M is required when argument PARTS is specified")') + case(583) + write (0,'("No more than one of the optional arguments ",a," must be present")')a_e_d case(600) write (0,'("Sparse Matrix and descriptors are in an invalid state for this subroutine call: ",i0)')i_e_d(1) case(700) diff --git a/base/newserial/psbn_base_mat_mod.f03 b/base/newserial/psbn_base_mat_mod.f03 index 72064c67..dc2bce07 100644 --- a/base/newserial/psbn_base_mat_mod.f03 +++ b/base/newserial/psbn_base_mat_mod.f03 @@ -51,6 +51,7 @@ module psbn_base_mat_mod procedure, pass(a) :: get_size procedure, pass(a) :: get_state procedure, pass(a) :: get_dupl + procedure, pass(a) :: get_fmt procedure, pass(a) :: is_null @@ -63,11 +64,10 @@ module psbn_base_mat_mod procedure, pass(a) :: is_triangle procedure, pass(a) :: is_unit procedure, pass(a) :: get_neigh - procedure, pass(a) :: allocate_mn procedure, pass(a) :: allocate_mnnz procedure, pass(a) :: reallocate_nz procedure, pass(a) :: free - generic, public :: allocate => allocate_mn, allocate_mnnz + generic, public :: allocate => allocate_mnnz generic, public :: reallocate => reallocate_nz procedure, pass(a) :: print => sparse_print @@ -80,10 +80,17 @@ 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 + & free, sparse_print,get_fmt contains + function get_fmt(a) result(res) + implicit none + class(psbn_base_sparse_mat), intent(in) :: a + character(len=5) :: res + res = 'NULL' + end function get_fmt + function get_dupl(a) result(res) implicit none class(psbn_base_sparse_mat), intent(in) :: a @@ -403,34 +410,12 @@ contains end subroutine get_neigh - subroutine allocate_mn(m,n,a) + subroutine allocate_mnnz(m,n,a,nz) use psb_error_mod implicit none integer, intent(in) :: m,n class(psbn_base_sparse_mat), intent(inout) :: a - - Integer :: err_act - character(len=20) :: name='allocate_mn' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - ! This is the base version. If we get here - ! it means the derived class is incomplete, - ! so we throw an error. - call psb_errpush(700,name) - - if (err_act /= psb_act_ret_) then - call psb_error() - end if - return - - end subroutine allocate_mn - - subroutine allocate_mnnz(m,n,nz,a) - use psb_error_mod - implicit none - integer, intent(in) :: m,n,nz - class(psbn_base_sparse_mat), intent(inout) :: a + integer, intent(in), optional :: nz Integer :: err_act character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 8e553d41..82cb7be5 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -43,20 +43,24 @@ module psbn_d_base_mat_mod procedure, pass(a) :: csins => d_coo_csins procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz procedure, pass(a) :: allocate_mnnz => d_coo_allocate_mnnz - procedure, pass(a) :: allocate_mn => d_coo_allocate_mn procedure, pass(a) :: cp_to_coo => d_cp_coo_to_coo procedure, pass(a) :: cp_from_coo => d_cp_coo_from_coo procedure, pass(a) :: cp_to_fmt => d_cp_coo_to_fmt procedure, pass(a) :: cp_from_fmt => d_cp_coo_from_fmt + procedure, pass(a) :: mv_to_coo => d_mv_coo_to_coo + procedure, pass(a) :: mv_from_coo => d_mv_coo_from_coo + procedure, pass(a) :: mv_to_fmt => d_mv_coo_to_fmt + procedure, pass(a) :: mv_from_fmt => d_mv_coo_from_fmt procedure, pass(a) :: fix => d_fix_coo procedure, pass(a) :: free => d_coo_free procedure, pass(a) :: print => d_coo_print + procedure, pass(a) :: get_fmt => d_coo_get_fmt end type psbn_d_coo_sparse_mat private :: d_coo_get_nzeros, d_coo_set_nzeros, & & d_coo_csmm, d_coo_csmv, d_coo_cssm, d_coo_cssv, & & d_coo_csins, d_coo_reallocate_nz, d_coo_allocate_mnnz, & - & d_coo_allocate_mn, d_fix_coo, d_coo_free, & + & 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 @@ -431,6 +435,12 @@ contains end subroutine mv_from_fmt + function d_coo_get_fmt(a) result(res) + implicit none + class(psbn_d_coo_sparse_mat), intent(in) :: a + character(len=5) :: res + res = 'COO' + end function d_coo_get_fmt subroutine d_fix_coo(a,info,idir) @@ -989,7 +999,6 @@ contains if (nz == 0) return nza = a%get_nzeros() -!!$ write(0,*) 'On entry to csins: ',nza call d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) if (info /= 0) goto 9999 @@ -1213,13 +1222,14 @@ contains end subroutine d_coo_free - subroutine d_coo_allocate_mnnz(m,n,nz,a) + subroutine d_coo_allocate_mnnz(m,n,a,nz) use psb_error_mod use psb_realloc_mod implicit none - integer, intent(in) :: m,n,nz + integer, intent(in) :: m,n class(psbn_d_coo_sparse_mat), intent(inout) :: a - Integer :: err_act, info + integer, intent(in), optional :: nz + Integer :: err_act, info, nz_ character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. @@ -1235,15 +1245,20 @@ contains call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) goto 9999 endif - if (nz < 0) then + if (present(nz)) then + nz_ = nz + else + nz_ = max(7*m,7*n,1) + end if + if (nz_ < 0) then info = 10 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) + if (info == 0) call psb_realloc(nz_,a%ia,info) + if (info == 0) call psb_realloc(nz_,a%ja,info) + if (info == 0) call psb_realloc(nz_,a%val,info) if (info == 0) then call a%set_nrows(m) call a%set_ncols(n) @@ -1267,47 +1282,6 @@ contains end subroutine d_coo_allocate_mnnz - - subroutine d_coo_allocate_mn(m,n,a) - use psb_error_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: m,n - class(psbn_d_coo_sparse_mat), intent(inout) :: a - Integer :: err_act, info, nz - character(len=20) :: name='allocate_mn' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - info = 0 - if (m < 0) then - info = 10 - call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) - goto 9999 - endif - if (n < 0) then - info = 10 - call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) - goto 9999 - endif - - nz = max(7*m,7*n,1) - call a%allocate(m,n,nz) - - 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_allocate_mn - subroutine d_coo_print(iout,a,iv,eirs,eics,head,ivr,ivc) use psb_spmat_type diff --git a/base/newserial/psbn_d_csr_impl.f03 b/base/newserial/psbn_d_csr_impl.f03 index 94836150..5c180c67 100644 --- a/base/newserial/psbn_d_csr_impl.f03 +++ b/base/newserial/psbn_d_csr_impl.f03 @@ -1226,7 +1226,6 @@ end subroutine d_cp_csr_from_coo_impl - subroutine d_cp_csr_to_coo_impl(a,b,info) use psb_const_mod use psbn_d_base_mat_mod @@ -1274,6 +1273,55 @@ subroutine d_cp_csr_to_coo_impl(a,b,info) end subroutine d_cp_csr_to_coo_impl +subroutine d_mv_csr_to_coo_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psbn_d_base_mat_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_to_coo_impl + implicit none + + class(psbn_d_csr_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + integer, allocatable :: itemp(:) + !locals + logical :: rwshr_ + Integer :: nza, nr, nc,i,j,irw, idl,err_act + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + nr = a%get_nrows() + 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 move_alloc(a%ja,b%ja) + call move_alloc(a%val,b%val) + call psb_realloc(nza,b%ia,info) + if (info /= 0) return + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + b%ia(j) = i + end do + end do + call a%free() + call b%fix(info) + + +end subroutine d_mv_csr_to_coo_impl + + subroutine d_mv_csr_from_coo_impl(a,b,info) use psb_const_mod @@ -1369,19 +1417,19 @@ subroutine d_mv_csr_from_coo_impl(a,b,info) end subroutine d_mv_csr_from_coo_impl -subroutine d_mv_csr_to_coo_impl(a,b,info) +subroutine d_mv_csr_to_fmt_impl(a,b,info) use psb_const_mod use psb_realloc_mod use psbn_d_base_mat_mod - use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_to_coo_impl + use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_to_fmt_impl implicit none class(psbn_d_csr_sparse_mat), intent(inout) :: a - class(psbn_d_coo_sparse_mat), intent(out) :: b + class(psbn_d_base_sparse_mat), intent(out) :: b integer, intent(out) :: info - integer, allocatable :: itemp(:) !locals + type(psbn_d_coo_sparse_mat) :: tmp logical :: rwshr_ Integer :: nza, nr, i,j,irw, idl,err_act, nc Integer, Parameter :: maxtry=8 @@ -1390,32 +1438,90 @@ subroutine d_mv_csr_to_coo_impl(a,b,info) info = 0 + call tmp%mv_from_fmt(a,info) + call b%mv_from_coo(tmp,info) - nr = a%get_nrows() - nc = a%get_ncols() - nza = 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()) - ! Dirty trick: call move_alloc to have the new data allocated just once. - call move_alloc(a%irp,itemp) - call move_alloc(a%ja,b%ja) - call move_alloc(a%val,b%val) - call psb_realloc(nza,b%ia,info) - if (info /= 0) return - if (allocated(itemp)) then - do i=1, nr - do j=itemp(i),itemp(i+1)-1 - b%ia(j) = i - end do - end do - end if - - call b%fix(info) +end subroutine d_mv_csr_to_fmt_impl -end subroutine d_mv_csr_to_coo_impl + +subroutine d_mv_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psbn_d_base_mat_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_mv_csr_from_fmt_impl + implicit none + + class(psbn_d_csr_sparse_mat), intent(inout) :: a + class(psbn_d_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + !locals + type(psbn_d_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + call tmp%mv_from_fmt(b,info) + call a%mv_from_coo(tmp,info) + +end subroutine d_mv_csr_from_fmt_impl + + + +subroutine d_cp_csr_to_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psbn_d_base_mat_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_cp_csr_to_fmt_impl + implicit none + + class(psbn_d_csr_sparse_mat), intent(in) :: a + class(psbn_d_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + !locals + type(psbn_d_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + call tmp%cp_from_fmt(a,info) + call b%mv_from_coo(tmp,info) + +end subroutine d_cp_csr_to_fmt_impl + + +subroutine d_cp_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psbn_d_base_mat_mod + use psbn_d_csr_mat_mod, psb_protect_name => d_cp_csr_from_fmt_impl + implicit none + + class(psbn_d_csr_sparse_mat), intent(inout) :: a + class(psbn_d_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + !locals + type(psbn_d_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + call tmp%cp_from_fmt(b,info) + call a%mv_from_coo(tmp,info) + +end subroutine d_cp_csr_from_fmt_impl diff --git a/base/newserial/psbn_d_csr_mat_mod.f03 b/base/newserial/psbn_d_csr_mat_mod.f03 index a1e493a5..f186d799 100644 --- a/base/newserial/psbn_d_csr_mat_mod.f03 +++ b/base/newserial/psbn_d_csr_mat_mod.f03 @@ -16,7 +16,6 @@ module psbn_d_csr_mat_mod procedure, pass(a) :: reallocate_nz => d_csr_reallocate_nz procedure, pass(a) :: csins => d_csr_csins procedure, pass(a) :: allocate_mnnz => d_csr_allocate_mnnz - procedure, pass(a) :: allocate_mn => d_csr_allocate_mn procedure, pass(a) :: cp_to_coo => d_cp_csr_to_coo procedure, pass(a) :: cp_from_coo => d_cp_csr_from_coo procedure, pass(a) :: cp_to_fmt => d_cp_csr_to_fmt @@ -27,12 +26,15 @@ module psbn_d_csr_mat_mod procedure, pass(a) :: mv_from_fmt => d_mv_csr_from_fmt procedure, pass(a) :: free => d_csr_free procedure, pass(a) :: print => d_csr_print + procedure, pass(a) :: get_fmt => d_csr_get_fmt 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_csins, d_csr_reallocate_nz, d_csr_allocate_mnnz, & - & d_csr_allocate_mn, d_csr_free, d_csr_print, & + & d_csr_free, d_csr_print, d_csr_get_fmt, & & d_cp_csr_to_coo, d_cp_csr_from_coo, & - & d_mv_csr_to_coo, d_mv_csr_from_coo + & d_mv_csr_to_coo, d_mv_csr_from_coo, & + & d_cp_csr_to_fmt, d_cp_csr_from_fmt, & + & d_mv_csr_to_fmt, d_mv_csr_from_fmt interface @@ -183,6 +185,14 @@ module psbn_d_csr_mat_mod contains + function d_csr_get_fmt(a) result(res) + implicit none + class(psbn_d_csr_sparse_mat), intent(in) :: a + character(len=5) :: res + res = 'CSR' + end function d_csr_get_fmt + + subroutine d_csr_reallocate_nz(nz,a) use psb_error_mod use psb_realloc_mod @@ -574,7 +584,7 @@ contains call psb_erractionsave(err_act) info = 0 -!!$ call d_cp_csr_to_fmt_impl(a,b,info) + call d_cp_csr_to_fmt_impl(a,b,info) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -606,7 +616,7 @@ contains call psb_erractionsave(err_act) info = 0 -!!$ call d_cp_csr_from_fmt_impl(a,b,info) + call d_cp_csr_from_fmt_impl(a,b,info) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -704,7 +714,7 @@ contains call psb_erractionsave(err_act) info = 0 -!!$ call d_mv_csr_to_fmt_impl(a,b,info) + call d_mv_csr_to_fmt_impl(a,b,info) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -736,7 +746,7 @@ contains call psb_erractionsave(err_act) info = 0 -!!$ call d_mv_csr_from_fmt_impl(a,b,info) + call d_mv_csr_from_fmt_impl(a,b,info) if (info /= 0) goto 9999 call psb_erractionrestore(err_act) @@ -755,13 +765,14 @@ contains end subroutine d_mv_csr_from_fmt - subroutine d_csr_allocate_mnnz(m,n,nz,a) + subroutine d_csr_allocate_mnnz(m,n,a,nz) use psb_error_mod use psb_realloc_mod implicit none - integer, intent(in) :: m,n,nz + integer, intent(in) :: m,n class(psbn_d_csr_sparse_mat), intent(inout) :: a - Integer :: err_act, info + integer, intent(in), optional :: nz + Integer :: err_act, info, nz_ character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. @@ -777,15 +788,20 @@ contains call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) goto 9999 endif - if (nz < 0) then + if (present(nz)) then + nz_ = nz + else + nz_ = max(7*m,7*n,1) + end if + if (nz_ < 0) then info = 10 call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) goto 9999 endif if (info == 0) call psb_realloc(m+1,a%irp,info) - if (info == 0) call psb_realloc(nz,a%ja,info) - if (info == 0) call psb_realloc(nz,a%val,info) + if (info == 0) call psb_realloc(nz_,a%ja,info) + if (info == 0) call psb_realloc(nz_,a%val,info) if (info == 0) then a%irp=0 call a%set_nrows(m) @@ -809,46 +825,6 @@ contains end subroutine d_csr_allocate_mnnz - subroutine d_csr_allocate_mn(m,n,a) - use psb_error_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: m,n - class(psbn_d_csr_sparse_mat), intent(inout) :: a - Integer :: err_act, info, nz - character(len=20) :: name='allocate_mn' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - info = 0 - if (m < 0) then - info = 10 - call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) - goto 9999 - endif - if (n < 0) then - info = 10 - call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) - goto 9999 - endif - - nz = max(7*m,7*n,1) - call a%allocate(m,n,nz) - - 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_allocate_mn - subroutine d_csr_print(iout,a,iv,eirs,eics,head,ivr,ivc) use psb_spmat_type use psb_string_mod diff --git a/base/newserial/psbn_mat_impl.f03 b/base/newserial/psbn_mat_impl.f03 index 9720e239..987c9ae7 100644 --- a/base/newserial/psbn_mat_impl.f03 +++ b/base/newserial/psbn_mat_impl.f03 @@ -1,3 +1,65 @@ +subroutine psbn_d_csall(nr,nc,a,info,nz) + use psbn_d_base_mat_mod + use psb_realloc_mod + use psb_sort_mod + use psbn_d_mat_mod, psb_protect_name => psbn_d_csall + implicit none + type(psbn_d_sparse_mat), intent(out) :: a + integer, intent(in) :: nr,nc + integer, intent(out) :: info + integer, intent(in), optional :: nz + + + info = 0 + call a%allocate(nr,nc,nz) + call a%set_state(psbn_spmat_bld_) + return + +end subroutine psbn_d_csall + + +subroutine psbn_d_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psbn_d_base_mat_mod + use psb_error_mod + use psbn_d_mat_mod, psb_protect_name => psbn_d_csins + implicit none + type(psbn_d_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + Integer :: err_act + character(len=20) :: name='psbn_csins' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.a%is_bld()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call a%a%csins(nz,val,ia,ja,imin,imax,jmin,jmax,info,gtl) + 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 psbn_d_csins + + + subroutine psbn_d_spcnv(a,b,info,type,mold,upd,dupl) use psbn_d_mat_mod, psb_protect_name => psbn_d_spcnv use psb_realloc_mod @@ -10,18 +72,98 @@ subroutine psbn_d_spcnv(a,b,info,type,mold,upd,dupl) character(len=*), optional, intent(in) :: type class(psbn_d_base_sparse_mat), intent(in), optional :: mold + + write(0,*) 'TO BE IMPLEMENTED ' + end subroutine psbn_d_spcnv subroutine psbn_d_spcnv_ip(a,info,type,mold,dupl) + use psb_error_mod + use psb_string_mod use psbn_d_mat_mod, psb_protect_name => psbn_d_spcnv_ip - use psb_realloc_mod - use psb_sort_mod implicit none - type(psbn_d_sparse_mat), intent(inout) :: a - integer, intent(out) :: info - integer,optional, intent(in) :: dupl - character(len=*), optional, intent(in) :: type + type(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer,optional, intent(in) :: dupl + character(len=*), optional, intent(in) :: type class(psbn_d_base_sparse_mat), intent(in), optional :: mold + + class(psbn_d_base_sparse_mat), allocatable :: altmp + class(psbn_d_base_sparse_mat), pointer :: aslct + type(psbn_d_csr_sparse_mat) :: csrtmp + Integer :: err_act + character(len=20) :: name='psbn_cscnv' + 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 + + if (present(dupl)) then + call a%set_dupl(dupl) + else + call a%set_dupl(psbn_dupl_def_) + end if + + if (count( (/present(mold),present(type) /)) > 1) then + info = 583 + call psb_errpush(info,name,a_err='TYPE, MOLD') + goto 9999 + end if + + if (present(mold)) then + + allocate(altmp, source=mold,stat=info) + + else if (present(type)) then + + select case (psb_toupper(type)) + case ('CSR') + allocate(psbn_d_csr_sparse_mat :: altmp, stat=info) + case ('COO') + allocate(psbn_d_coo_sparse_mat :: altmp, stat=info) + case default + info = 136 + call psb_errpush(info,name,a_err=type) + goto 9999 + end select + else + allocate(altmp, source=csrtmp,stat=info) + end if + + select type ( aa => a%a ) + class is (psbn_d_coo_sparse_mat) + ! Quick route from coo + call altmp%mv_from_coo(aa, info) + class default + call altmp%mv_from_fmt(aa, info) + end select + + if (info /= 0) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call move_alloc(altmp,a%a) + call a%set_asb() + + 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 psbn_d_spcnv_ip diff --git a/base/newserial/psbn_mat_mod.f03 b/base/newserial/psbn_mat_mod.f03 index c72f91ce..32449b3a 100644 --- a/base/newserial/psbn_mat_mod.f03 +++ b/base/newserial/psbn_mat_mod.f03 @@ -2,13 +2,28 @@ module psbn_d_mat_mod use psbn_d_base_mat_mod - + use psbn_d_csr_mat_mod + type :: psbn_d_sparse_mat class(psbn_d_base_sparse_mat), allocatable :: a contains + procedure, pass(a) :: set_nrows + procedure, pass(a) :: set_ncols + procedure, pass(a) :: set_dupl + procedure, pass(a) :: set_state + procedure, pass(a) :: set_null + procedure, pass(a) :: set_bld + procedure, pass(a) :: set_upd + procedure, pass(a) :: set_asb + procedure, pass(a) :: set_sorted + procedure, pass(a) :: set_upper + procedure, pass(a) :: set_lower + procedure, pass(a) :: set_triangle + procedure, pass(a) :: set_unit + procedure, pass(a) :: get_nrows procedure, pass(a) :: get_ncols procedure, pass(a) :: get_nzeros @@ -26,11 +41,13 @@ module psbn_d_mat_mod procedure, pass(a) :: is_triangle procedure, pass(a) :: is_unit procedure, pass(a) :: get_neigh - procedure, pass(a) :: allocate_mn procedure, pass(a) :: allocate_mnnz procedure, pass(a) :: reallocate_nz procedure, pass(a) :: free - generic, public :: allocate => allocate_mn, allocate_mnnz + procedure, pass(a) :: print => sparse_print + procedure, pass(a) :: get_fmt => sparse_get_fmt + + generic, public :: allocate => allocate_mnnz generic, public :: reallocate => reallocate_nz procedure, pass(a) :: d_csmv @@ -46,18 +63,45 @@ module psbn_d_mat_mod private :: get_nrows, get_ncols, get_nzeros, get_size, & & get_state, get_dupl, is_null, is_bld, is_upd, & & is_asb, is_sorted, is_upper, is_lower, is_triangle, & - & is_unit, get_neigh, allocate_mn, allocate_mnnz, & - & reallocate_nz, free, d_csmv, d_csmm, d_cssv, d_cssm + & is_unit, get_neigh, allocate_mnnz, & + & reallocate_nz, free, d_csmv, d_csmm, d_cssv, d_cssm, sparse_print, & + & set_nrows, set_ncols, set_dupl, set_state, set_null, set_bld, & + & set_upd, set_asb, set_sorted, set_upper, set_lower, set_triangle, & + & set_unit + - interface psbn_spcnv + interface psbn_csall + subroutine psbn_d_csall(nr,nc,a,info,nz) + use psbn_d_base_mat_mod + import psbn_d_sparse_mat + type(psbn_d_sparse_mat), intent(out) :: a + integer, intent(in) :: nr,nc + integer, intent(out) :: info + integer, intent(in), optional :: nz + end subroutine psbn_d_csall + end interface + + interface psbn_csins + subroutine psbn_d_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psbn_d_base_mat_mod + import psbn_d_sparse_mat + type(psbn_d_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + end subroutine psbn_d_csins + end interface + + interface psbn_cscnv subroutine psbn_d_spcnv(a,b,info,type,mold,upd,dupl) use psbn_d_base_mat_mod import psbn_d_sparse_mat type(psbn_d_sparse_mat), intent(in) :: a type(psbn_d_sparse_mat), intent(out) :: b integer, intent(out) :: info - integer,optional, intent(in) :: dupl, upd + integer, optional, intent(in) :: dupl, upd character(len=*), optional, intent(in) :: type class(psbn_d_base_sparse_mat), intent(in), optional :: mold @@ -75,6 +119,22 @@ module psbn_d_mat_mod contains + + function sparse_get_fmt(a) result(res) + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + character(len=5) :: res + + if (allocated(a%a)) then + res = a%a%get_fmt() + else + res = 'NULL' + end if + + end function sparse_get_fmt + + + function get_dupl(a) result(res) use psb_error_mod implicit none @@ -244,13 +304,44 @@ contains end function is_sorted + + subroutine set_nrows(m,a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(in) :: m + Integer :: err_act, info + character(len=20) :: name='set_nrows' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif - function get_nzeros(a) result(res) + call a%a%set_nrows(m) + + 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 set_nrows + + subroutine set_ncols(n,a) use psb_error_mod implicit none - class(psbn_d_sparse_mat), intent(in) :: a - integer :: res - + class(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(in) :: n Integer :: err_act, info character(len=20) :: name='get_nzeros' logical, parameter :: debug=.false. @@ -261,8 +352,39 @@ contains call psb_errpush(info,name) goto 9999 endif + call a%a%set_ncols(n) - res = a%a%get_nzeros() + 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 set_ncols + + + subroutine set_state(n,a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + call a%a%set_state(n) call psb_erractionrestore(err_act) return @@ -275,16 +397,17 @@ contains return end if - end function get_nzeros - function get_size(a) result(res) + end subroutine set_state + + + subroutine set_dupl(n,a) use psb_error_mod implicit none - class(psbn_d_sparse_mat), intent(in) :: a - integer :: res - + class(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(in) :: n Integer :: err_act, info - character(len=20) :: name='get_size' + character(len=20) :: name='get_nzeros' logical, parameter :: debug=.false. call psb_erractionsave(err_act) @@ -293,8 +416,8 @@ contains call psb_errpush(info,name) goto 9999 endif - - res = a%a%get_size() + + call a%a%set_dupl(n) call psb_erractionrestore(err_act) return @@ -306,26 +429,79 @@ contains call psb_error() return end if + + + end subroutine set_dupl + + subroutine set_null(a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_null() + + call psb_erractionrestore(err_act) return - end function get_size +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if - subroutine get_neigh(a,idx,neigh,n,info,lev) + end subroutine set_null + + subroutine set_bld(a) use psb_error_mod implicit none - class(psbn_d_sparse_mat), intent(in) :: a - integer, intent(in) :: idx - integer, intent(out) :: n - integer, allocatable, intent(out) :: neigh(:) - integer, intent(out) :: info - integer, optional, intent(in) :: lev + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_bld() - Integer :: err_act - character(len=20) :: name='get_neigh' + 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 set_bld + + subroutine set_upd(a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' logical, parameter :: debug=.false. - info = 0 call psb_erractionsave(err_act) if (.not.allocated(a%a)) then info = 1121 @@ -333,9 +509,38 @@ contains goto 9999 endif - call a%a%get_neigh(idx,neigh,n,info,lev) + call a%a%set_upd() - 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 set_upd + + subroutine set_asb(a) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_asb() call psb_erractionrestore(err_act) return @@ -347,64 +552,218 @@ contains call psb_error() return end if + + end subroutine set_asb + + subroutine set_sorted(a,val) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_sorted(val) + + call psb_erractionrestore(err_act) return - end subroutine get_neigh +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + end subroutine set_sorted - subroutine allocate_mn(m,n,a,type,mold) + subroutine set_triangle(a,val) use psb_error_mod - use psb_string_mod implicit none - integer, intent(in) :: m,n class(psbn_d_sparse_mat), intent(inout) :: a - character(len=*), intent(in), optional :: type - class(psbn_d_base_sparse_mat), intent(in), optional :: mold + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_triangle(val) + + 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 set_triangle + + subroutine set_unit(a,val) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val Integer :: err_act, info - character(len=20) :: name='allocate_mn' - character(len=8) :: type_ + character(len=20) :: name='get_nzeros' logical, parameter :: debug=.false. + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_unit(val) + + 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 set_unit + + subroutine set_lower(a,val) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. call psb_erractionsave(err_act) - info = 0 - if (allocated(a%a)) then - call a%a%free() - deallocate(a%a) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_lower(val) + + 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 - if (present(mold)) then - allocate(a%a, source=mold, stat=info) + end subroutine set_lower - else + subroutine set_upper(a,val) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. - if (present(type)) then - type_ = psb_toupper(type) - else - type_ = 'COO' - end if + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif - select case(type) - case('COO') - allocate(psbn_d_coo_sparse_mat :: a%a, stat=info) -! Add here a few other data structures inplemented by default. + call a%a%set_upper(val) -!!$ case('CSR') -!!$ allocate(psbn_d_csr_sparse_mat :: a%a, stat=info) + call psb_erractionrestore(err_act) + return - case default - allocate(psbn_d_coo_sparse_mat :: a%a, stat=info) - end select +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return end if - if (info /= 0) then - info = 4010 + end subroutine set_upper + + + + + function get_nzeros(a) result(res) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) goto 9999 + endif + + res = a%a%get_nzeros() + + 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 function get_nzeros + + function get_size(a) result(res) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res - call a%a%allocate(m,n) - + Integer :: err_act, info + character(len=20) :: name='get_size' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%get_size() + call psb_erractionrestore(err_act) return @@ -417,15 +776,96 @@ contains end if return + end function get_size - end subroutine allocate_mn + subroutine sparse_print(iout,a,iv,eirs,eics,head,ivr,ivc) + use psb_error_mod + implicit none + + integer, intent(in) :: iout + class(psbn_d_sparse_mat), intent(in) :: a + integer, intent(in), optional :: iv(:) + integer, intent(in), optional :: eirs,eics + character(len=*), optional :: head + integer, intent(in), optional :: ivr(:), ivc(:) + + Integer :: err_act, info + character(len=20) :: name='sparse_print' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif - subroutine allocate_mnnz(m,n,nz,a,type,mold) + call a%a%print(iout,iv,eirs,eics,head,ivr,ivc) + + + 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 sparse_print + + subroutine get_neigh(a,idx,neigh,n,info,lev) + use psb_error_mod + implicit none + class(psbn_d_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer, intent(out) :: n + integer, allocatable, intent(out) :: neigh(:) + integer, intent(out) :: info + integer, optional, intent(in) :: lev + + Integer :: err_act + character(len=20) :: name='get_neigh' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%get_neigh(idx,neigh,n,info,lev) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine get_neigh + + + subroutine allocate_mnnz(m,n,a,nz,type,mold) use psb_error_mod use psb_string_mod implicit none - integer, intent(in) :: m,n,nz + integer, intent(in) :: m,n class(psbn_d_sparse_mat), intent(inout) :: a + integer, intent(in), optional :: nz character(len=*), intent(in), optional :: type class(psbn_d_base_sparse_mat), intent(in), optional :: mold diff --git a/test/serial/Makefile b/test/serial/Makefile index 5c949230..10f0ff42 100644 --- a/test/serial/Makefile +++ b/test/serial/Makefile @@ -14,10 +14,13 @@ FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). EXEDIR=./runs -all: d_coo_matgen +all: d_coo_matgen d_matgen d_coo_matgen: d_coo_matgen.o $(F90LINK) d_coo_matgen.o -o d_coo_matgen $(PSBLAS_LIB) $(LDLIBS) /bin/mv d_coo_matgen $(EXEDIR) +d_matgen: d_matgen.o + $(F90LINK) d_matgen.o -o d_matgen $(PSBLAS_LIB) $(LDLIBS) + /bin/mv d_matgen $(EXEDIR) #ppde spde @@ -34,7 +37,7 @@ spde: spde.o clean: - /bin/rm -f d_coo_matgen.o tpg.o ppde.o spde.o $(EXEDIR)/ppde + /bin/rm -f d_coo_matgen.o d_matgen.o tpg.o ppde.o spde.o $(EXEDIR)/ppde verycleanlib: (cd ../..; make veryclean) lib: diff --git a/test/serial/d_matgen.f03 b/test/serial/d_matgen.f03 new file mode 100644 index 00000000..2a5c13aa --- /dev/null +++ b/test/serial/d_matgen.f03 @@ -0,0 +1,469 @@ +! +program d_coo_matgen + use psb_base_mod + use psb_prec_mod + use psb_krylov_mod + use psbn_d_base_mat_mod + use psbn_d_csr_mat_mod + use psbn_d_mat_mod + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer :: idim + + ! miscellaneous + real(psb_dpk_), parameter :: one = 1.d0 + real(psb_dpk_) :: t1, t2, tprec + + ! sparse matrix and preconditioner + type(psb_dspmat_type) :: a + type(psb_dprec_type) :: prec + ! descriptor + type(psb_desc_type) :: desc_a + ! dense matrices + real(psb_dpk_), allocatable :: b(:), x(:) + ! blacs parameters + integer :: ictxt, iam, np + + ! solver parameters + integer :: iter, itmax,itrace, istopc, irst + integer(psb_long_int_k_) :: amatsize, precsize, descsize + real(psb_dpk_) :: err, eps + + ! other variables + integer :: info, err_act + character(len=20) :: name,ch_err + + info=0 + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + if(psb_get_errstatus() /= 0) goto 9999 + + call psb_set_errverbosity(2) + + ! + ! get parameters + ! + call get_parms(ictxt,idim) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= 0) then + call psb_error(ictxt) + end if + + call psb_exit(ictxt) + stop + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + end if + stop + +contains + ! + ! get iteration parameters from standard input + ! + subroutine get_parms(ictxt,idim) + integer :: ictxt + integer :: idim + integer :: np, iam + integer :: intbuf(10), ip + + call psb_info(ictxt, iam, np) + + read(*,*) idim + + + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: pde90 methd prec dim & + &[istop itmax itrace]' + write(iout,*)' where:' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' prec : bjac diag none' + write(iout,*)' dim number of points along each axis' + write(iout,*)' the size of the resulting linear ' + write(iout,*)' system is dim**3' + write(iout,*)' istop stopping criterion 1, 2 ' + write(iout,*)' itmax maximum number of iterations [500] ' + write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' >= 1 do tracing every itrace' + write(iout,*)' iterations ' + end subroutine pr_usage + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,afmt,info) + ! + ! discretize the partial diferential equation + ! + ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) + ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. + ! + ! Boundary conditions are set in a very simple way, by adding + ! equations of the form + ! + ! u(x,y) = exp(-x^2-y^2-z^2) + ! + ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. + ! + use psb_base_mod + implicit none + integer :: idim + integer, parameter :: nb=20 + real(psb_dpk_), allocatable :: b(:),xv(:) + type(psb_desc_type) :: desc_a + integer :: ictxt, info + character :: afmt*5 + type(psb_dspmat_type) :: a + real(psb_dpk_) :: zt(nb),glob_x,glob_y,glob_z + integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer :: x,y,z,ia,indx_owner + integer :: np, iam, nr, nt,nz,isz + integer :: element + integer, allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + type(psbn_d_sparse_mat) :: a_n + type(psbn_d_coo_sparse_mat) :: acoo + type(psbn_d_csr_sparse_mat) :: acsr + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah + real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcpy, tmov + real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 + external :: a1, a2, a3, a4, b1, b2, b3 + integer :: err_act + + character(len=20) :: name, ch_err + + info = 0 + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + deltah = 1.d0/(idim-1) + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(0,*) iam, 'Initialization error ',nr,nt,m + write(0,*) iam, 'Initialization ',nr,nt,m + nlr = nt + call psb_barrier(ictxt) + + t0 = psb_wtime() + + call psbn_csall(nr,nr,a_n,info) +!!$ call acoo%allocate + + talc = psb_wtime()-t0 + +!!$ write(*,*) 'Test get size:',d_coo_get_size(acoo) +!!$ write(*,*) 'Test 2 get size:',acoo%get_size(),acoo%get_nzeros() + + if (info /= 0) then + info=4010 + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= 0 ) then + info=4000 + call psb_errpush(info,name) + goto 9999 + endif + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) +!!$ write(0,*) 'Row ',ii,ib + element = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=i + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + x = glob_row/(idim*idim) + else + x = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then + y = (glob_row-(x-1)*idim*idim)/idim + else + y = (glob_row-(x-1)*idim*idim)/idim+1 + endif + z = glob_row-(x-1)*idim*idim-(y-1)*idim + ! glob_x, glob_y, glob_x coordinates + glob_x=x*deltah + glob_y=y*deltah + glob_z=z*deltah + + ! check on boundary points + zt(k) = 0.d0 + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + if (x==1) then + val(element)=-b1(glob_x,glob_y,glob_z)& + & -a1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-glob_y**2-glob_z**2)*(-val(element)) + else + val(element)=-b1(glob_x,glob_y,glob_z)& + & -a1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-2)*idim*idim+(y-1)*idim+(z) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y-1,z) + if (y==1) then + val(element)=-b2(glob_x,glob_y,glob_z)& + & -a2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b2(glob_x,glob_y,glob_z)& + & -a2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-1)*idim*idim+(y-2)*idim+(z) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z-1) + if (z==1) then + val(element)=-b3(glob_x,glob_y,glob_z)& + & -a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b3(glob_x,glob_y,glob_z)& + & -a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-1)*idim*idim+(y-1)*idim+(z-1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z) + val(element)=2*b1(glob_x,glob_y,glob_z)& + & +2*b2(glob_x,glob_y,glob_z)& + & +2*b3(glob_x,glob_y,glob_z)& + & +a1(glob_x,glob_y,glob_z)& + & +a2(glob_x,glob_y,glob_z)& + & +a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-1)*idim*idim+(y-1)*idim+(z) + irow(element) = glob_row + element = element+1 + ! term depending on (x,y,z+1) + if (z==idim) then + val(element)=-b1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-1)*idim*idim+(y-1)*idim+(z+1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y+1,z) + if (y==idim) then + val(element)=-b2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (x-1)*idim*idim+(y)*idim+(z) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x+1,y,z) + if (x