fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-18 11:52:11 +00:00
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
Changelog
 base/modules/psb_base_mat_mod.f03
 base/modules/psb_base_mod.f90
 base/modules/psb_d_base_mat_mod.f03
 base/modules/psb_d_csr_mat_mod.f03
 base/modules/psb_inter_desc_mod.f90
 base/modules/psb_inter_desc_type.f90
 base/modules/psb_linmap_mod.f90
 base/modules/psb_linmap_type_mod.f90
 base/modules/psb_mat_mod.f03
 base/modules/psb_spmat_type.f03
 base/modules/psb_tools_mod.f90
 base/psblas/psb_dnrmi.f90
 base/psblas/psb_dspmm.f90
 base/psblas/psb_dspsm.f90
 base/serial/Makefile
 base/serial/coo/Makefile
 base/serial/csr/Makefile
 base/serial/dp/Makefile
 base/serial/f03/psbn_d_coo_impl.f03
 base/serial/f03/psbn_d_csr_impl.f03
 base/serial/f77/Makefile
 base/serial/jad/Makefile
 base/serial/psb_getrow_mod.f90
 base/serial/psb_regen_mod.f90
 base/serial/psb_update_mod.f90
 base/tools/psb_dcdbldext.F90
 base/tools/psb_dspalloc.f90
 base/tools/psb_dspasb.f90
 base/tools/psb_dspfree.f90
 base/tools/psb_dsphalo.F90
 base/tools/psb_dspins.f90
 base/tools/psb_dsprn.f90
 base/tools/psb_linmap.f90
 krylov/psb_dbicg.f90
 krylov/psb_dcg.F90
 krylov/psb_dcgs.f90
 krylov/psb_dcgstab.F90
 krylov/psb_dcgstabl.f90
 krylov/psb_drgmres.f90
 krylov/psb_krylov_mod.f90
 prec/psb_dbjac_aply.f90
 prec/psb_dbjac_bld.f90
 prec/psb_ddiagsc_bld.f90
 prec/psb_dilu_fct.f90
 prec/psb_dprecbld.f90
 prec/psb_prec_mod.f90
 test/fileread/df_sample.f90
 test/fileread/runs/dfs.inp
 test/pargen/ppde.f90
 util/psb_hbio_mod.f90
 util/psb_mat_dist_mod.f90
 util/psb_mmio_mod.f90

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