mirror of
https://github.com/sfilippone/psblas3.git
synced 2026-10-24 14:51:49 +00:00
Merged MixedI8 in new branch (to be later merged into development)
This commit is contained in:
Salvatore Filippone
@@ -36,14 +36,14 @@ module psb_blockpart_mod
|
||||
|
||||
contains
|
||||
subroutine part_block(global_indx,n,np,pv,nv)
|
||||
use psb_base_mod, only : psb_ipk_, psb_mpik_
|
||||
use psb_base_mod, only : psb_ipk_, psb_mpk_, psb_lpk_
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), intent(in) :: global_indx, n
|
||||
integer(psb_lpk_), intent(in) :: global_indx, n
|
||||
integer(psb_ipk_), intent(in) :: np
|
||||
integer(psb_ipk_), intent(out) :: nv
|
||||
integer(psb_ipk_), intent(out) :: pv(*)
|
||||
integer(psb_ipk_) :: dim_block
|
||||
integer(psb_lpk_) :: dim_block
|
||||
|
||||
dim_block = (n + np - 1)/np
|
||||
nv = 1
|
||||
@@ -56,10 +56,11 @@ contains
|
||||
|
||||
|
||||
subroutine bld_partblock(n,np,ivg)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
integer(psb_ipk_) :: n,np,ivg(*)
|
||||
use psb_base_mod, only : psb_ipk_, psb_mpk_, psb_lpk_
|
||||
integer(psb_lpk_) :: n
|
||||
integer(psb_ipk_) :: np,ivg(*)
|
||||
|
||||
integer(psb_ipk_) :: dim_block,i
|
||||
integer(psb_lpk_) :: dim_block,i
|
||||
|
||||
|
||||
dim_block = (n + np - 1)/np
|
||||
@@ -69,7 +70,5 @@ contains
|
||||
|
||||
end subroutine bld_partblock
|
||||
|
||||
|
||||
|
||||
end module psb_blockpart_mod
|
||||
|
||||
|
||||
@@ -361,3 +361,338 @@ subroutine chb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine chb_write
|
||||
|
||||
|
||||
subroutine lchb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
complex(psb_spk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
|
||||
character :: rhstype*3,type*3,key*8
|
||||
character(len=72) :: mtitle_
|
||||
character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
|
||||
integer(psb_lpk_) :: indcrd, ptrcrd, totcrd,&
|
||||
& valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
|
||||
type(psb_lc_csc_sparse_mat) :: acsc
|
||||
type(psb_lc_coo_sparse_mat) :: acoo
|
||||
integer(psb_ipk_) :: ircode, infile, info
|
||||
integer(psb_lpk_) :: i,nzr
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
iret = 0
|
||||
ircode = 0
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
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 acsc%allocate(nrow,ncol,nnzero)
|
||||
if (ircode /= 0 ) then
|
||||
write(psb_err_unit,*) 'Memory allocation failed'
|
||||
goto 993
|
||||
end if
|
||||
|
||||
if (present(mtitle)) mtitle=mtitle_
|
||||
|
||||
|
||||
if (psb_tolower(type(1:1)) == 'c') then
|
||||
if (psb_tolower(type(2:2)) == 'u') then
|
||||
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
call a%mv_from(acsc)
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 's') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 'h') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = conjg(acoo%val(i))
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
|
||||
call a%cscnv(ircode,type='csr')
|
||||
if (infile /= 5) close(infile)
|
||||
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 iret=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 iret=902
|
||||
write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
|
||||
return
|
||||
993 iret=993
|
||||
write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
|
||||
return
|
||||
end subroutine lchb_read
|
||||
|
||||
subroutine lchb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(in), target :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
complex(psb_spk_), optional :: rhs(:), g(:), x(:)
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
character(len=*), parameter:: ptrfmt='(10I8)',indfmt='(10I8)'
|
||||
integer(psb_ipk_), parameter :: jptr=10,jind=10
|
||||
character(len=*), parameter:: valfmt='(4E20.12)',rhsfmt='(4E20.12)'
|
||||
integer(psb_ipk_), parameter :: jval=2,jrhs=2
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
type(psb_lc_csc_sparse_mat) :: acsc
|
||||
character(len=72) :: mtitle_
|
||||
character(len=8) :: key_
|
||||
|
||||
character :: rhstype*3,type*3
|
||||
|
||||
integer(psb_lpk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
|
||||
& nrow,ncol,nnzero, neltvl, nrhs, nrhsix
|
||||
|
||||
iret = 0
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
if (present(mtitle)) then
|
||||
mtitle_ = mtitle
|
||||
else
|
||||
mtitle_ = 'Temporary PSBLAS title '
|
||||
endif
|
||||
if (present(key)) then
|
||||
key_ = key
|
||||
else
|
||||
key_ = 'PSBMAT00'
|
||||
endif
|
||||
|
||||
|
||||
call acsc%cp_from_fmt(a%a, iret)
|
||||
if (iret /= 0) return
|
||||
|
||||
|
||||
nrow = acsc%get_nrows()
|
||||
ncol = acsc%get_ncols()
|
||||
nnzero = acsc%get_nzeros()
|
||||
|
||||
neltvl = 0
|
||||
|
||||
ptrcrd = (ncol+1)/jptr
|
||||
if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
|
||||
indcrd = nnzero/jind
|
||||
if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
|
||||
valcrd = nnzero/jval
|
||||
if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
|
||||
rhstype = ''
|
||||
if (present(rhs)) then
|
||||
if (size(rhs)<nrow) then
|
||||
rhscrd = 0
|
||||
else
|
||||
rhscrd = nrow/jrhs
|
||||
if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
|
||||
endif
|
||||
nrhs = 1
|
||||
rhstype(1:1) = 'F'
|
||||
else
|
||||
rhscrd = 0
|
||||
nrhs = 0
|
||||
end if
|
||||
totcrd = ptrcrd + indcrd + valcrd + rhscrd
|
||||
|
||||
nrhsix = nrhs*nrow
|
||||
|
||||
if (present(g)) then
|
||||
rhstype(2:2) = 'G'
|
||||
end if
|
||||
if (present(x)) then
|
||||
rhstype(3:3) = 'X'
|
||||
end if
|
||||
type = 'CUA'
|
||||
|
||||
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) (acsc%icp(i),i=1,ncol+1)
|
||||
write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) write (iout,fmt=valfmt) (acsc%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)
|
||||
|
||||
|
||||
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
iret=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine lchb_write
|
||||
|
||||
@@ -88,11 +88,11 @@ subroutine psb_cmatdist(a_glob, a, ictxt, desc_a,&
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing
|
||||
integer(psb_ipk_) :: i_count, j_count,&
|
||||
& k_count, root, liwork, nrow, ncol, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act, int_err(5)
|
||||
integer(psb_ipk_) :: k_count, root, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_ipk_), allocatable :: irow(:),icol(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
complex(psb_spk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
@@ -140,8 +140,7 @@ subroutine psb_cmatdist(a_glob, a, ictxt, desc_a,&
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
int_err(1)=liwork
|
||||
call psb_errpush(info,name,i_err=int_err,a_err='integer')
|
||||
call psb_errpush(info,name,i_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
@@ -339,3 +338,319 @@ subroutine psb_cmatdist(a_glob, a, ictxt, desc_a,&
|
||||
return
|
||||
|
||||
end subroutine psb_cmatdist
|
||||
|
||||
|
||||
|
||||
subroutine psb_lcmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_cspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_cspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
use psb_base_mod
|
||||
use psb_mat_mod
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lcspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_cspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_c_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing, root, iproc
|
||||
integer(psb_ipk_) :: err_act, il, inz
|
||||
integer(psb_lpk_) :: k_count, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, nnr, err
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
complex(psb_spk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info = psb_success_
|
||||
err = 0
|
||||
name = 'psb_c_mat_dist'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
! executable statements
|
||||
if (present(inroot)) then
|
||||
root = inroot
|
||||
else
|
||||
root = psb_root_
|
||||
end if
|
||||
call psb_info(ictxt, iam, np)
|
||||
if (iam == root) then
|
||||
nrow = a_glob%get_nrows()
|
||||
ncol = a_glob%get_ncols()
|
||||
if (nrow /= ncol) then
|
||||
write(psb_err_unit,*) 'a rectangular matrix ? ',nrow,ncol
|
||||
info=-1
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
endif
|
||||
nnzero = a_glob%get_nzeros()
|
||||
nrhs = 1
|
||||
endif
|
||||
|
||||
use_parts = present(parts)
|
||||
use_v = present(v)
|
||||
if (count((/ use_parts, use_v /)) /= 1) then
|
||||
info=psb_err_no_optional_arg_
|
||||
call psb_errpush(info,name,a_err=" v, parts")
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
! broadcast informations to other processors
|
||||
call psb_bcast(ictxt,nrow, root)
|
||||
call psb_bcast(ictxt,ncol, root)
|
||||
call psb_bcast(ictxt,nnzero, root)
|
||||
call psb_bcast(ictxt,nrhs, root)
|
||||
liwork = max(np, nrow + ncol)
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,l_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
write (*, fmt = *) 'start matdist',root, size(iwork),&
|
||||
&nrow, ncol, nnzero,nrhs
|
||||
endif
|
||||
if (use_parts) then
|
||||
call psb_cdall(ictxt,desc_a,info,mg=nrow,parts=parts)
|
||||
else if (use_v) then
|
||||
call psb_cdall(ictxt,desc_a,info,vg=v)
|
||||
else
|
||||
info = -1
|
||||
end if
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
inz = ((nnzero+np-1)/np)
|
||||
call psb_spall(a,desc_a,info,nnz=inz)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
isize = 3*nb*max(((nnzero+nrow)/nrow),nb)
|
||||
allocate(val(isize),irow(isize),icol(isize),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
i_count = 1
|
||||
|
||||
do while (i_count <= nrow)
|
||||
|
||||
if (use_parts) then
|
||||
call parts(i_count,nrow,np,iwork, np_sharing)
|
||||
!
|
||||
! np_sharing allows for overlap in the data distribution.
|
||||
! If an index is overlapped, then we have to send its row
|
||||
! to multiple processes. NOTE: we are assuming the output
|
||||
! from PARTS is coherent, otherwise a deadlock is the most
|
||||
! likely outcome.
|
||||
!
|
||||
j_count = i_count
|
||||
if (np_sharing == 1) then
|
||||
iproc = iwork(1)
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
call parts(j_count,nrow,np,iwrk2, np_sharing)
|
||||
if (np_sharing /= 1 ) exit
|
||||
if (iwrk2(1) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
else
|
||||
np_sharing = 1
|
||||
j_count = i_count
|
||||
iproc = v(i_count)
|
||||
iwork(1:np_sharing) = iproc
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
if (v(j_count) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
|
||||
! now we should insert rows i_count..j_count-1
|
||||
nnr = j_count - i_count
|
||||
|
||||
if (iam == root) then
|
||||
|
||||
ll = 0
|
||||
do i= i_count, j_count-1
|
||||
call a_glob%csget(i,i,nz,&
|
||||
& irow,icol,val,info,nzin=ll,append=.true.)
|
||||
if (info /= psb_success_) then
|
||||
if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then
|
||||
write(psb_err_unit,*) 'Allocation failure? This should not happen!'
|
||||
end if
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
ll = ll + nz
|
||||
end do
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
|
||||
if (iproc == iam) then
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
else
|
||||
call psb_snd(ictxt,nnr,iproc)
|
||||
call psb_snd(ictxt,ll,iproc)
|
||||
call psb_snd(ictxt,irow(1:ll),iproc)
|
||||
call psb_snd(ictxt,icol(1:ll),iproc)
|
||||
call psb_snd(ictxt,val(1:ll),iproc)
|
||||
call psb_rcv(ictxt,ll,iproc)
|
||||
endif
|
||||
end do
|
||||
else if (iam /= root) then
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
if (iproc == iam) then
|
||||
call psb_rcv(ictxt,nnr,root)
|
||||
call psb_rcv(ictxt,ll,root)
|
||||
if (ll > size(irow)) then
|
||||
write(psb_err_unit,*) iam,'need to reallocate ',ll
|
||||
deallocate(val,irow,icol)
|
||||
allocate(val(ll),irow(ll),icol(ll),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
endif
|
||||
call psb_rcv(ictxt,irow(1:ll),root)
|
||||
call psb_rcv(ictxt,icol(1:ll),root)
|
||||
call psb_rcv(ictxt,val(1:ll),root)
|
||||
call psb_snd(ictxt,ll,root)
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psspins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
end do
|
||||
endif
|
||||
i_count = j_count
|
||||
end do
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t0 = psb_wtime()
|
||||
call psb_cdasb(desc_a,info)
|
||||
t1 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime()
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=fmt,mold=mold)
|
||||
t3 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (iam == root) then
|
||||
write(psb_out_unit,*) 'descriptor assembly: ',t1-t0
|
||||
write(psb_out_unit,*) 'sparse matrix assembly: ',t3-t2
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(val,irow,icol,iwork,iwrk2,stat=info)
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='deallocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (iam == root) write (*, fmt = *) 'end matdist'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ictxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
end subroutine psb_lcmatdist
|
||||
|
||||
@@ -31,7 +31,8 @@
|
||||
!
|
||||
module psb_c_mat_dist_mod
|
||||
use psb_base_mod, only : psb_ipk_, psb_spk_, psb_desc_type, psb_parts, &
|
||||
& psb_cspmat_type, psb_c_base_sparse_mat, psb_c_vect_type
|
||||
& psb_cspmat_type, psb_c_base_sparse_mat, psb_c_vect_type, &
|
||||
& psb_lcspmat_type
|
||||
|
||||
interface psb_matdist
|
||||
subroutine psb_cmatdist(a_glob, a, ictxt, desc_a,&
|
||||
@@ -91,6 +92,64 @@ module psb_c_mat_dist_mod
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_cmatdist
|
||||
subroutine psb_lcmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_lcspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_cspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
import :: psb_ipk_, psb_cspmat_type, psb_spk_, psb_desc_type,&
|
||||
& psb_c_base_sparse_mat, psb_c_vect_type, psb_parts, &
|
||||
& psb_lcspmat_type
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lcspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_cspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_c_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_lcmatdist
|
||||
end interface
|
||||
|
||||
end module psb_c_mat_dist_mod
|
||||
|
||||
@@ -461,3 +461,179 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine cmm_mat_write
|
||||
|
||||
subroutine lcmm_mat_read(a, info, iunit, filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character :: mmheader*15, fmt*15, object*10, type*10, sym*15
|
||||
character(1024) :: line
|
||||
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr
|
||||
integer(psb_ipk_) :: ircode,infile
|
||||
type(psb_lc_coo_sparse_mat), allocatable :: acoo
|
||||
real(psb_spk_) :: are, aim
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
read(infile,fmt=*,end=902) mmheader, object, fmt, type, sym
|
||||
|
||||
if ( (psb_tolower(object) /= 'matrix').or.(psb_tolower(fmt) /= 'coordinate')) then
|
||||
write(psb_err_unit,*) 'READ_MATRIX: input file type not yet supported'
|
||||
info=909
|
||||
return
|
||||
end if
|
||||
|
||||
do
|
||||
read(infile,fmt='(a)') line
|
||||
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) == 'complex').and.(psb_tolower(sym) == 'general')) then
|
||||
call acoo%allocate(nrow,ncol,nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
|
||||
end do
|
||||
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) == 'complex').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 acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(info)
|
||||
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
|
||||
else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
call acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_spk_)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = conjg(acoo%val(i))
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(info)
|
||||
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
info=904
|
||||
end if
|
||||
if (infile /= 5) close(infile)
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 info=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 info=902
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Unexpected end of file '
|
||||
return
|
||||
905 info=905
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Error at line',i
|
||||
return
|
||||
993 info=993
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Memory allocation failure'
|
||||
return
|
||||
end subroutine lcmm_mat_read
|
||||
|
||||
|
||||
subroutine lcmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
call a%print(iout,head=mtitle)
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
info=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine lcmm_mat_write
|
||||
|
||||
@@ -267,7 +267,7 @@ contains
|
||||
name = 'mat_renum_amd'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
#if defined(HAVE_AMD)
|
||||
#if defined(HAVE_AMD) && defined(IPK4)
|
||||
|
||||
info = psb_success_
|
||||
nr = a%get_nrows()
|
||||
|
||||
@@ -314,3 +314,290 @@ subroutine dhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine dhb_write
|
||||
|
||||
subroutine ldhb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
real(psb_dpk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
|
||||
character :: rhstype*3,type*3,key*8
|
||||
character(len=72) :: mtitle_
|
||||
character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
|
||||
integer(psb_lpk_) :: indcrd, ptrcrd, totcrd,&
|
||||
& valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
|
||||
type(psb_ld_csc_sparse_mat) :: acsc
|
||||
type(psb_ld_coo_sparse_mat) :: acoo
|
||||
integer(psb_ipk_) :: ircode, infile, info
|
||||
integer(psb_lpk_) :: i,nzr
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
iret = 0
|
||||
ircode = 0
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
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 acsc%allocate(nrow,ncol,nnzero)
|
||||
if (ircode /= 0 ) then
|
||||
write(psb_err_unit,*) 'Memory allocation failed'
|
||||
goto 993
|
||||
end if
|
||||
|
||||
if (present(mtitle)) mtitle=mtitle_
|
||||
|
||||
|
||||
if (psb_tolower(type(1:1)) == 'r') then
|
||||
if (psb_tolower(type(2:2)) == 'u') then
|
||||
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
call a%mv_from(acsc)
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 's') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
|
||||
call a%cscnv(ircode,type='csr')
|
||||
if (infile /= 5) close(infile)
|
||||
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 iret=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 iret=902
|
||||
write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
|
||||
return
|
||||
993 iret=993
|
||||
write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
|
||||
return
|
||||
end subroutine ldhb_read
|
||||
|
||||
subroutine ldhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(in), target :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
real(psb_dpk_), optional :: rhs(:), g(:), x(:)
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
character(len=*), parameter:: ptrfmt='(10I8)',indfmt='(10I8)'
|
||||
integer(psb_ipk_), parameter :: jptr=10,jind=10
|
||||
character(len=*), parameter:: valfmt='(4E20.12)',rhsfmt='(4E20.12)'
|
||||
integer(psb_ipk_), parameter :: jval=4,jrhs=4
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
type(psb_ld_csc_sparse_mat) :: acsc
|
||||
character(len=72) :: mtitle_
|
||||
character(len=8) :: key_
|
||||
|
||||
character :: rhstype*3,type*3
|
||||
|
||||
integer(psb_lpk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
|
||||
& nrow,ncol,nnzero, neltvl, nrhs, nrhsix
|
||||
|
||||
iret = 0
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
if (present(mtitle)) then
|
||||
mtitle_ = mtitle
|
||||
else
|
||||
mtitle_ = 'Temporary PSBLAS title '
|
||||
endif
|
||||
if (present(key)) then
|
||||
key_ = key
|
||||
else
|
||||
key_ = 'PSBMAT00'
|
||||
endif
|
||||
|
||||
|
||||
call acsc%cp_from_fmt(a%a, iret)
|
||||
if (iret /= 0) return
|
||||
|
||||
|
||||
nrow = acsc%get_nrows()
|
||||
ncol = acsc%get_ncols()
|
||||
nnzero = acsc%get_nzeros()
|
||||
|
||||
neltvl = 0
|
||||
|
||||
ptrcrd = (ncol+1)/jptr
|
||||
if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
|
||||
indcrd = nnzero/jind
|
||||
if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
|
||||
valcrd = nnzero/jval
|
||||
if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
|
||||
rhstype = ''
|
||||
if (present(rhs)) then
|
||||
if (size(rhs)<nrow) then
|
||||
rhscrd = 0
|
||||
else
|
||||
rhscrd = nrow/jrhs
|
||||
if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
|
||||
endif
|
||||
nrhs = 1
|
||||
rhstype(1:1) = 'F'
|
||||
else
|
||||
rhscrd = 0
|
||||
nrhs = 0
|
||||
end if
|
||||
totcrd = ptrcrd + indcrd + valcrd + rhscrd
|
||||
|
||||
nrhsix = nrhs*nrow
|
||||
|
||||
if (present(g)) then
|
||||
rhstype(2:2) = 'G'
|
||||
end if
|
||||
if (present(x)) then
|
||||
rhstype(3:3) = 'X'
|
||||
end if
|
||||
type = 'RUA'
|
||||
|
||||
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) (acsc%icp(i),i=1,ncol+1)
|
||||
write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) write (iout,fmt=valfmt) (acsc%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)
|
||||
|
||||
|
||||
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
iret=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine ldhb_write
|
||||
|
||||
@@ -88,11 +88,11 @@ subroutine psb_dmatdist(a_glob, a, ictxt, desc_a,&
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing
|
||||
integer(psb_ipk_) :: i_count, j_count,&
|
||||
& k_count, root, liwork, nrow, ncol, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act, int_err(5)
|
||||
integer(psb_ipk_) :: k_count, root, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_ipk_), allocatable :: irow(:),icol(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
real(psb_dpk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
@@ -140,8 +140,7 @@ subroutine psb_dmatdist(a_glob, a, ictxt, desc_a,&
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
int_err(1)=liwork
|
||||
call psb_errpush(info,name,i_err=int_err,a_err='integer')
|
||||
call psb_errpush(info,name,i_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
@@ -339,3 +338,319 @@ subroutine psb_dmatdist(a_glob, a, ictxt, desc_a,&
|
||||
return
|
||||
|
||||
end subroutine psb_dmatdist
|
||||
|
||||
|
||||
|
||||
subroutine psb_ldmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_dspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_dspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
use psb_base_mod
|
||||
use psb_mat_mod
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_ldspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_dspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_d_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing, root, iproc
|
||||
integer(psb_ipk_) :: err_act, il, inz
|
||||
integer(psb_lpk_) :: k_count, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, nnr, err
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
real(psb_dpk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info = psb_success_
|
||||
err = 0
|
||||
name = 'psb_d_mat_dist'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
! executable statements
|
||||
if (present(inroot)) then
|
||||
root = inroot
|
||||
else
|
||||
root = psb_root_
|
||||
end if
|
||||
call psb_info(ictxt, iam, np)
|
||||
if (iam == root) then
|
||||
nrow = a_glob%get_nrows()
|
||||
ncol = a_glob%get_ncols()
|
||||
if (nrow /= ncol) then
|
||||
write(psb_err_unit,*) 'a rectangular matrix ? ',nrow,ncol
|
||||
info=-1
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
endif
|
||||
nnzero = a_glob%get_nzeros()
|
||||
nrhs = 1
|
||||
endif
|
||||
|
||||
use_parts = present(parts)
|
||||
use_v = present(v)
|
||||
if (count((/ use_parts, use_v /)) /= 1) then
|
||||
info=psb_err_no_optional_arg_
|
||||
call psb_errpush(info,name,a_err=" v, parts")
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
! broadcast informations to other processors
|
||||
call psb_bcast(ictxt,nrow, root)
|
||||
call psb_bcast(ictxt,ncol, root)
|
||||
call psb_bcast(ictxt,nnzero, root)
|
||||
call psb_bcast(ictxt,nrhs, root)
|
||||
liwork = max(np, nrow + ncol)
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,l_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
write (*, fmt = *) 'start matdist',root, size(iwork),&
|
||||
&nrow, ncol, nnzero,nrhs
|
||||
endif
|
||||
if (use_parts) then
|
||||
call psb_cdall(ictxt,desc_a,info,mg=nrow,parts=parts)
|
||||
else if (use_v) then
|
||||
call psb_cdall(ictxt,desc_a,info,vg=v)
|
||||
else
|
||||
info = -1
|
||||
end if
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
inz = ((nnzero+np-1)/np)
|
||||
call psb_spall(a,desc_a,info,nnz=inz)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
isize = 3*nb*max(((nnzero+nrow)/nrow),nb)
|
||||
allocate(val(isize),irow(isize),icol(isize),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
i_count = 1
|
||||
|
||||
do while (i_count <= nrow)
|
||||
|
||||
if (use_parts) then
|
||||
call parts(i_count,nrow,np,iwork, np_sharing)
|
||||
!
|
||||
! np_sharing allows for overlap in the data distribution.
|
||||
! If an index is overlapped, then we have to send its row
|
||||
! to multiple processes. NOTE: we are assuming the output
|
||||
! from PARTS is coherent, otherwise a deadlock is the most
|
||||
! likely outcome.
|
||||
!
|
||||
j_count = i_count
|
||||
if (np_sharing == 1) then
|
||||
iproc = iwork(1)
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
call parts(j_count,nrow,np,iwrk2, np_sharing)
|
||||
if (np_sharing /= 1 ) exit
|
||||
if (iwrk2(1) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
else
|
||||
np_sharing = 1
|
||||
j_count = i_count
|
||||
iproc = v(i_count)
|
||||
iwork(1:np_sharing) = iproc
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
if (v(j_count) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
|
||||
! now we should insert rows i_count..j_count-1
|
||||
nnr = j_count - i_count
|
||||
|
||||
if (iam == root) then
|
||||
|
||||
ll = 0
|
||||
do i= i_count, j_count-1
|
||||
call a_glob%csget(i,i,nz,&
|
||||
& irow,icol,val,info,nzin=ll,append=.true.)
|
||||
if (info /= psb_success_) then
|
||||
if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then
|
||||
write(psb_err_unit,*) 'Allocation failure? This should not happen!'
|
||||
end if
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
ll = ll + nz
|
||||
end do
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
|
||||
if (iproc == iam) then
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
else
|
||||
call psb_snd(ictxt,nnr,iproc)
|
||||
call psb_snd(ictxt,ll,iproc)
|
||||
call psb_snd(ictxt,irow(1:ll),iproc)
|
||||
call psb_snd(ictxt,icol(1:ll),iproc)
|
||||
call psb_snd(ictxt,val(1:ll),iproc)
|
||||
call psb_rcv(ictxt,ll,iproc)
|
||||
endif
|
||||
end do
|
||||
else if (iam /= root) then
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
if (iproc == iam) then
|
||||
call psb_rcv(ictxt,nnr,root)
|
||||
call psb_rcv(ictxt,ll,root)
|
||||
if (ll > size(irow)) then
|
||||
write(psb_err_unit,*) iam,'need to reallocate ',ll
|
||||
deallocate(val,irow,icol)
|
||||
allocate(val(ll),irow(ll),icol(ll),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
endif
|
||||
call psb_rcv(ictxt,irow(1:ll),root)
|
||||
call psb_rcv(ictxt,icol(1:ll),root)
|
||||
call psb_rcv(ictxt,val(1:ll),root)
|
||||
call psb_snd(ictxt,ll,root)
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psspins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
end do
|
||||
endif
|
||||
i_count = j_count
|
||||
end do
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t0 = psb_wtime()
|
||||
call psb_cdasb(desc_a,info)
|
||||
t1 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime()
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=fmt,mold=mold)
|
||||
t3 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (iam == root) then
|
||||
write(psb_out_unit,*) 'descriptor assembly: ',t1-t0
|
||||
write(psb_out_unit,*) 'sparse matrix assembly: ',t3-t2
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(val,irow,icol,iwork,iwrk2,stat=info)
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='deallocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (iam == root) write (*, fmt = *) 'end matdist'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ictxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
end subroutine psb_ldmatdist
|
||||
|
||||
@@ -31,7 +31,8 @@
|
||||
!
|
||||
module psb_d_mat_dist_mod
|
||||
use psb_base_mod, only : psb_ipk_, psb_dpk_, psb_desc_type, psb_parts, &
|
||||
& psb_dspmat_type, psb_d_base_sparse_mat, psb_d_vect_type
|
||||
& psb_dspmat_type, psb_d_base_sparse_mat, psb_d_vect_type, &
|
||||
& psb_ldspmat_type
|
||||
|
||||
interface psb_matdist
|
||||
subroutine psb_dmatdist(a_glob, a, ictxt, desc_a,&
|
||||
@@ -91,6 +92,64 @@ module psb_d_mat_dist_mod
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_dmatdist
|
||||
subroutine psb_ldmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_ldspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_dspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
import :: psb_ipk_, psb_dspmat_type, psb_dpk_, psb_desc_type,&
|
||||
& psb_d_base_sparse_mat, psb_d_vect_type, psb_parts, &
|
||||
& psb_ldspmat_type
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_ldspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_dspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_d_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_ldmatdist
|
||||
end interface
|
||||
|
||||
end module psb_d_mat_dist_mod
|
||||
|
||||
@@ -452,3 +452,180 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine dmm_mat_write
|
||||
|
||||
|
||||
subroutine ldmm_mat_read(a, info, iunit, filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character :: mmheader*15, fmt*15, object*10, type*10, sym*15
|
||||
character(1024) :: line
|
||||
integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr
|
||||
integer(psb_ipk_) :: ircode, infile
|
||||
type(psb_ld_coo_sparse_mat), allocatable :: acoo
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
read(infile,fmt=*,end=902) mmheader, object, fmt, type, sym
|
||||
|
||||
if ( (psb_tolower(object) /= 'matrix').or.(psb_tolower(fmt) /= 'coordinate')) then
|
||||
write(psb_err_unit,*) 'READ_MATRIX: input file type not yet supported'
|
||||
info=909
|
||||
return
|
||||
end if
|
||||
|
||||
do
|
||||
read(infile,fmt='(a)') line
|
||||
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 acoo%allocate(nrow,ncol,nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
|
||||
end do
|
||||
call acoo%set_nzeros(nnzero)
|
||||
|
||||
else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then
|
||||
call acoo%allocate(nrow,ncol,nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
|
||||
end do
|
||||
acoo%val(:) = done
|
||||
call acoo%set_nzeros(nnzero)
|
||||
|
||||
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 acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
|
||||
else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then
|
||||
call acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i)
|
||||
end do
|
||||
acoo%val(:) = done
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
info=904
|
||||
end if
|
||||
|
||||
if (info == 0) then
|
||||
call acoo%fix(info)
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
end if
|
||||
|
||||
if (infile /= 5) close(infile)
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 info=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 info=902
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Unexpected end of file '
|
||||
return
|
||||
905 info=905
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Error at line',i
|
||||
return
|
||||
993 info=993
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Memory allocation failure'
|
||||
return
|
||||
end subroutine ldmm_mat_read
|
||||
|
||||
|
||||
subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
call a%print(iout,head=mtitle)
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
info=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine ldmm_mat_write
|
||||
|
||||
|
||||
|
||||
@@ -267,7 +267,7 @@ contains
|
||||
name = 'mat_renum_amd'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
#if defined(HAVE_AMD)
|
||||
#if defined(HAVE_AMD) && defined(IPK4)
|
||||
|
||||
info = psb_success_
|
||||
nr = a%get_nrows()
|
||||
|
||||
+88
-1
@@ -33,7 +33,10 @@ module psb_hbio_mod
|
||||
|
||||
use psb_base_mod, only : psb_ipk_, psb_spk_, psb_dpk_,&
|
||||
& psb_sspmat_type, psb_cspmat_type, &
|
||||
& psb_dspmat_type, psb_zspmat_type
|
||||
& psb_dspmat_type, psb_zspmat_type, &
|
||||
& psb_lsspmat_type, psb_lcspmat_type, &
|
||||
& psb_ldspmat_type, psb_lzspmat_type
|
||||
|
||||
|
||||
|
||||
public hb_read, hb_write
|
||||
@@ -78,6 +81,46 @@ module psb_hbio_mod
|
||||
complex(psb_dpk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
end subroutine zhb_read
|
||||
subroutine lshb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
import :: psb_lsspmat_type, psb_spk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
real(psb_spk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
end subroutine lshb_read
|
||||
subroutine ldhb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
import :: psb_ldspmat_type, psb_dpk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
real(psb_dpk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
end subroutine ldhb_read
|
||||
subroutine lchb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
import :: psb_lcspmat_type, psb_spk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
complex(psb_spk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
end subroutine lchb_read
|
||||
subroutine lzhb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
import :: psb_lzspmat_type, psb_dpk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
complex(psb_dpk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
end subroutine lzhb_read
|
||||
end interface
|
||||
|
||||
interface hb_write
|
||||
@@ -125,6 +168,50 @@ module psb_hbio_mod
|
||||
character(len=*), optional, intent(in) :: key
|
||||
complex(psb_dpk_), optional :: rhs(:), g(:), x(:)
|
||||
end subroutine zhb_write
|
||||
subroutine lshb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
import :: psb_lsspmat_type, psb_spk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(inout) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
real(psb_spk_), optional :: rhs(:), g(:), x(:)
|
||||
end subroutine lshb_write
|
||||
subroutine ldhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
import :: psb_ldspmat_type, psb_dpk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(inout) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
real(psb_dpk_), optional :: rhs(:), g(:), x(:)
|
||||
end subroutine ldhb_write
|
||||
subroutine lchb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
import :: psb_lcspmat_type, psb_spk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(inout) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
complex(psb_spk_), optional :: rhs(:), g(:), x(:)
|
||||
end subroutine lchb_write
|
||||
subroutine lzhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
import :: psb_lzspmat_type, psb_dpk_, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(inout) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
complex(psb_dpk_), optional :: rhs(:), g(:), x(:)
|
||||
end subroutine lzhb_write
|
||||
end interface
|
||||
|
||||
end module psb_hbio_mod
|
||||
|
||||
@@ -54,8 +54,9 @@
|
||||
! uses information prepared by the previous two subroutines.
|
||||
!
|
||||
module psb_metispart_mod
|
||||
use psb_base_mod, only : psb_ipk_, psb_sspmat_type, psb_cspmat_type,&
|
||||
& psb_dspmat_type, psb_zspmat_type, psb_err_unit, psb_mpik_,&
|
||||
use psb_base_mod, only : psb_sspmat_type, psb_cspmat_type,&
|
||||
& psb_dspmat_type, psb_zspmat_type, psb_err_unit, &
|
||||
& psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_, &
|
||||
& psb_s_csr_sparse_mat, psb_d_csr_sparse_mat, &
|
||||
& psb_c_csr_sparse_mat, psb_z_csr_sparse_mat
|
||||
public part_graph, build_mtpart, distr_mtpart,&
|
||||
@@ -76,7 +77,7 @@ contains
|
||||
|
||||
subroutine part_graph(global_indx,n,np,pv,nv)
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(in) :: global_indx, n
|
||||
integer(psb_lpk_), intent(in) :: global_indx, n
|
||||
integer(psb_ipk_), intent(in) :: np
|
||||
integer(psb_ipk_), intent(out) :: nv
|
||||
integer(psb_ipk_), intent(out) :: pv(*)
|
||||
@@ -101,8 +102,9 @@ contains
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: root, ictxt
|
||||
integer(psb_ipk_) :: n, me, np, info
|
||||
|
||||
integer(psb_ipk_) :: me, np, info
|
||||
integer(psb_lpk_) :: n
|
||||
|
||||
call psb_info(ictxt,me,np)
|
||||
|
||||
if (.not.((root>=0).and.(root<np))) then
|
||||
@@ -325,7 +327,7 @@ contains
|
||||
integer(psb_ipk_), allocatable :: irpl(:),jal(:),gvl(:)
|
||||
real(psb_spk_),allocatable :: wgh_(:)
|
||||
|
||||
#if defined(HAVE_METIS)
|
||||
#if defined(HAVE_METIS) && defined(IPK4)
|
||||
interface
|
||||
! subroutine METIS_PartGraphKway(n,ixadj,iadj,ivwg,iajw,&
|
||||
! & wgflag,numflag,nparts,weights,iopt,nedc,part) bind(c)
|
||||
|
||||
+71
-1
@@ -34,7 +34,9 @@ module psb_mmio_mod
|
||||
|
||||
use psb_base_mod, only : psb_ipk_, psb_spk_, psb_dpk_,&
|
||||
& psb_sspmat_type, psb_cspmat_type, &
|
||||
& psb_dspmat_type, psb_zspmat_type
|
||||
& psb_dspmat_type, psb_zspmat_type, &
|
||||
& psb_lsspmat_type, psb_lcspmat_type, &
|
||||
& psb_ldspmat_type, psb_lzspmat_type
|
||||
|
||||
public mm_mat_read, mm_mat_write, mm_array_read, mm_array_write
|
||||
|
||||
@@ -270,6 +272,38 @@ module psb_mmio_mod
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine zmm_mat_read
|
||||
subroutine lsmm_mat_read(a, info, iunit, filename)
|
||||
import :: psb_lsspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lsmm_mat_read
|
||||
subroutine ldmm_mat_read(a, info, iunit, filename)
|
||||
import :: psb_ldspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine ldmm_mat_read
|
||||
subroutine lcmm_mat_read(a, info, iunit, filename)
|
||||
import :: psb_lcspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lcmm_mat_read
|
||||
subroutine lzmm_mat_read(a, info, iunit, filename)
|
||||
import :: psb_lzspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lzmm_mat_read
|
||||
end interface
|
||||
|
||||
interface mm_mat_write
|
||||
@@ -309,6 +343,42 @@ module psb_mmio_mod
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine zmm_mat_write
|
||||
subroutine lsmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
import :: psb_lsspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lsmm_mat_write
|
||||
subroutine ldmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
import :: psb_ldspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_ldspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine ldmm_mat_write
|
||||
subroutine lcmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
import :: psb_lcspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lcspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lcmm_mat_write
|
||||
subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
import :: psb_lzspmat_type, psb_ipk_
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
end subroutine lzmm_mat_write
|
||||
end interface
|
||||
|
||||
end module psb_mmio_mod
|
||||
|
||||
+359
-30
@@ -38,7 +38,7 @@
|
||||
! See test/pargen/psb_X_pdeNd for examples of usage
|
||||
!
|
||||
module psb_partidx_mod
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
|
||||
interface idx2ijk
|
||||
module procedure idx2ijk3d, idx2ijkv, idx2ijk2d
|
||||
@@ -47,6 +47,15 @@ module psb_partidx_mod
|
||||
interface ijk2idx
|
||||
module procedure ijk2idx3d, ijk2idxv, ijk2idx2d
|
||||
end interface ijk2idx
|
||||
interface idx2ijk
|
||||
module procedure lidx2ijk3d, lidx2ijkv, lidx2ijk2d,&
|
||||
& lidx2lijk3d, lidx2lijkv, lidx2lijk2d
|
||||
end interface idx2ijk
|
||||
|
||||
interface ijk2idx
|
||||
module procedure ijk2lidx3d, ijk2lidxv, ijk2lidx2d,&
|
||||
& lijk2lidx3d, lijk2lidxv, lijk2lidx2d
|
||||
end interface ijk2idx
|
||||
|
||||
logical, private, save :: col_major=.true.
|
||||
|
||||
@@ -85,13 +94,13 @@ contains
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine idx2ijk3d(i,j,k,idx,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: i,j,k
|
||||
integer(psb_ipk_), intent(in) :: idx,nx,ny,nz
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(out) :: i,j,k
|
||||
integer(psb_mpk_), intent(in) :: idx,nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_ipk_) :: coords(3)
|
||||
integer(psb_mpk_) :: coords(3)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny,nz],base)
|
||||
|
||||
@@ -102,13 +111,13 @@ contains
|
||||
end subroutine idx2ijk3d
|
||||
|
||||
subroutine idx2ijk2d(i,j,idx,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: i,j
|
||||
integer(psb_ipk_), intent(in) :: idx,nx,ny
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(out) :: i,j
|
||||
integer(psb_mpk_), intent(in) :: idx,nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_ipk_) :: coords(2)
|
||||
integer(psb_mpk_) :: coords(2)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny],base)
|
||||
|
||||
@@ -130,13 +139,13 @@ contains
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine idx2ijkv(coords,idx,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: coords(:)
|
||||
integer(psb_ipk_), intent(in) :: idx,dims(:)
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(out) :: coords(:)
|
||||
integer(psb_mpk_), intent(in) :: idx,dims(:)
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_ipk_) :: base_, idx_, i, sz
|
||||
integer(psb_mpk_) :: base_, idx_, i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
@@ -170,7 +179,188 @@ contains
|
||||
end if
|
||||
|
||||
end subroutine idx2ijkv
|
||||
|
||||
subroutine lidx2ijk3d(i,j,k,idx,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_mpk_), intent(out) :: i,j,k
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_mpk_), intent(in) :: nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_mpk_) :: coords(3)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny,nz],base)
|
||||
|
||||
k = coords(3)
|
||||
j = coords(2)
|
||||
i = coords(1)
|
||||
|
||||
end subroutine lidx2ijk3d
|
||||
|
||||
subroutine lidx2ijk2d(i,j,idx,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_mpk_), intent(out) :: i,j
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_mpk_), intent(in) :: nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_mpk_) :: coords(2)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny],base)
|
||||
|
||||
j = coords(2)
|
||||
i = coords(1)
|
||||
|
||||
end subroutine lidx2ijk2d
|
||||
|
||||
!
|
||||
! Given a global index IDX and the domain size (NX,NY,NZ)
|
||||
! compute the point coordinates (I,J,K)
|
||||
! Optional argument: base 0 or 1, default 1
|
||||
!
|
||||
! This mapping is equivalent to a loop nesting:
|
||||
! idx = base
|
||||
! do i=1,nx
|
||||
! do j=1,ny
|
||||
! do k=1,nz
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine lidx2ijkv(coords,idx,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_mpk_), intent(out) :: coords(:)
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_mpk_), intent(in) :: dims(:)
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_epk_) :: base_, idx_
|
||||
integer(psb_mpk_) :: i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
base_ = 1
|
||||
end if
|
||||
|
||||
idx_ = idx - base_
|
||||
|
||||
if (size(coords) < size(dims)) then
|
||||
write(0,*) 'Error: size mismatch ',size(coords),size(dims)
|
||||
coords = 0
|
||||
return
|
||||
end if
|
||||
|
||||
!
|
||||
! This code is equivalent to (3D case)
|
||||
! k = mod(idx_,nz) + base_
|
||||
! j = mod(idx_/nz,ny) + base_
|
||||
! i = mod(idx_/(nx*ny),nx) + base_
|
||||
!
|
||||
if (col_major) then
|
||||
do i=1,size(dims)
|
||||
coords(i) = mod(idx_,dims(i)) + base_
|
||||
idx_ = idx_ / dims(i)
|
||||
end do
|
||||
else
|
||||
do i=size(dims),1,-1
|
||||
coords(i) = mod(idx_,dims(i)) + base_
|
||||
idx_ = idx_ / dims(i)
|
||||
end do
|
||||
end if
|
||||
|
||||
end subroutine lidx2ijkv
|
||||
|
||||
subroutine lidx2lijk3d(i,j,k,idx,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: i,j,k
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_epk_), intent(in) :: nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_epk_) :: coords(3)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny,nz],base)
|
||||
|
||||
k = coords(3)
|
||||
j = coords(2)
|
||||
i = coords(1)
|
||||
|
||||
end subroutine lidx2lijk3d
|
||||
|
||||
subroutine lidx2lijk2d(i,j,idx,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: i,j
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_epk_), intent(in) :: nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_epk_) :: coords(2)
|
||||
|
||||
call idx2ijk(coords,idx,[nx,ny],base)
|
||||
|
||||
j = coords(2)
|
||||
i = coords(1)
|
||||
|
||||
end subroutine lidx2lijk2d
|
||||
|
||||
!
|
||||
! Given a global index IDX and the domain size (NX,NY,NZ)
|
||||
! compute the point coordinates (I,J,K)
|
||||
! Optional argument: base 0 or 1, default 1
|
||||
!
|
||||
! This mapping is equivalent to a loop nesting:
|
||||
! idx = base
|
||||
! do i=1,nx
|
||||
! do j=1,ny
|
||||
! do k=1,nz
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine lidx2lijkv(coords,idx,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: coords(:)
|
||||
integer(psb_epk_), intent(in) :: idx
|
||||
integer(psb_epk_), intent(in) :: dims(:)
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_epk_) :: base_, idx_
|
||||
integer(psb_epk_) :: i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
base_ = 1
|
||||
end if
|
||||
|
||||
idx_ = idx - base_
|
||||
|
||||
if (size(coords) < size(dims)) then
|
||||
write(0,*) 'Error: size mismatch ',size(coords),size(dims)
|
||||
coords = 0
|
||||
return
|
||||
end if
|
||||
|
||||
!
|
||||
! This code is equivalent to (3D case)
|
||||
! k = mod(idx_,nz) + base_
|
||||
! j = mod(idx_/nz,ny) + base_
|
||||
! i = mod(idx_/(nx*ny),nx) + base_
|
||||
!
|
||||
if (col_major) then
|
||||
do i=1,size(dims)
|
||||
coords(i) = mod(idx_,dims(i)) + base_
|
||||
idx_ = idx_ / dims(i)
|
||||
end do
|
||||
else
|
||||
do i=size(dims),1,-1
|
||||
coords(i) = mod(idx_,dims(i)) + base_
|
||||
idx_ = idx_ / dims(i)
|
||||
end do
|
||||
end if
|
||||
|
||||
end subroutine lidx2lijkv
|
||||
!
|
||||
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
|
||||
! compute the global index IDX
|
||||
@@ -184,13 +374,13 @@ contains
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine ijk2idxv(idx,coords,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(in) :: coords(:),dims(:)
|
||||
integer(psb_ipk_), intent(out) :: idx
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(in) :: coords(:),dims(:)
|
||||
integer(psb_mpk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_ipk_) :: base_, i, sz
|
||||
integer(psb_mpk_) :: base_, i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
@@ -230,27 +420,166 @@ contains
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine ijk2idx3d(idx,i,j,k,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: idx
|
||||
integer(psb_ipk_), intent(in) :: i,j,k,nx,ny,nz
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in) :: i,j,k,nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*nz*ny + (j-base_)*nz + k - base_) + base_
|
||||
call ijk2idx(idx,[i,j,k],[nx,ny,nz],base)
|
||||
end subroutine ijk2idx3d
|
||||
|
||||
subroutine ijk2idx2d(idx,i,j,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: idx
|
||||
integer(psb_ipk_), intent(in) :: i,j,nx,ny
|
||||
integer(psb_ipk_), intent(in), optional :: base
|
||||
integer(psb_mpk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in) :: i,j,nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*ny + (j-base_) + base_
|
||||
call ijk2idx(idx,[i,j],[nx,ny],base)
|
||||
end subroutine ijk2idx2d
|
||||
|
||||
subroutine ijk2lidxv(idx,coords,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_mpk_), intent(in) :: coords(:),dims(:)
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_mpk_) :: base_, i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
base_ = 1
|
||||
end if
|
||||
sz = size(coords)
|
||||
if (sz /= size(dims)) then
|
||||
write(0,*) 'Error: size mismatch ',size(coords),size(dims)
|
||||
idx = 0
|
||||
return
|
||||
end if
|
||||
|
||||
if (col_major) then
|
||||
idx = coords(sz) - base_
|
||||
do i=sz-1,1,-1
|
||||
idx = (idx * dims(i)) + coords(i) - base_
|
||||
end do
|
||||
idx = idx + base_
|
||||
else
|
||||
idx = coords(1) - base_
|
||||
do i=2,sz
|
||||
idx = (idx * dims(i)) + coords(i) - base_
|
||||
end do
|
||||
idx = idx + base_
|
||||
end if
|
||||
|
||||
end subroutine ijk2lidxv
|
||||
!
|
||||
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
|
||||
! compute the global index IDX
|
||||
! Optional argument: base 0 or 1, default 1
|
||||
!
|
||||
! This mapping is equivalent to a loop nesting:
|
||||
! idx = base
|
||||
! do i=1,nx
|
||||
! do j=1,ny
|
||||
! do k=1,nz
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine ijk2lidx3d(idx,i,j,k,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in) :: i,j,k,nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*nz*ny + (j-base_)*nz + k - base_) + base_
|
||||
call ijk2idx(idx,[i,j,k],[nx,ny,nz],base)
|
||||
end subroutine ijk2lidx3d
|
||||
|
||||
subroutine ijk2lidx2d(idx,i,j,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in) :: i,j,nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*ny + (j-base_) + base_
|
||||
call ijk2idx(idx,[i,j],[nx,ny],base)
|
||||
end subroutine ijk2lidx2d
|
||||
|
||||
|
||||
subroutine lijk2lidxv(idx,coords,dims,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(in) :: coords(:),dims(:)
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
integer(psb_epk_) :: base_, i, sz
|
||||
if (present(base)) then
|
||||
base_ = base
|
||||
else
|
||||
base_ = 1
|
||||
end if
|
||||
sz = size(coords)
|
||||
if (sz /= size(dims)) then
|
||||
write(0,*) 'Error: size mismatch ',size(coords),size(dims)
|
||||
idx = 0
|
||||
return
|
||||
end if
|
||||
|
||||
if (col_major) then
|
||||
idx = coords(sz) - base_
|
||||
do i=sz-1,1,-1
|
||||
idx = (idx * dims(i)) + coords(i) - base_
|
||||
end do
|
||||
idx = idx + base_
|
||||
else
|
||||
idx = coords(1) - base_
|
||||
do i=2,sz
|
||||
idx = (idx * dims(i)) + coords(i) - base_
|
||||
end do
|
||||
idx = idx + base_
|
||||
end if
|
||||
|
||||
end subroutine lijk2lidxv
|
||||
!
|
||||
! Given a triple (I,J,K) and the domain size (NX,NY,NZ)
|
||||
! compute the global index IDX
|
||||
! Optional argument: base 0 or 1, default 1
|
||||
!
|
||||
! This mapping is equivalent to a loop nesting:
|
||||
! idx = base
|
||||
! do i=1,nx
|
||||
! do j=1,ny
|
||||
! do k=1,nz
|
||||
! ijk2idx(i,j,k) = idx
|
||||
! idx = idx + 1
|
||||
subroutine lijk2lidx3d(idx,i,j,k,nx,ny,nz,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_epk_), intent(in) :: i,j,k,nx,ny,nz
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*nz*ny + (j-base_)*nz + k - base_) + base_
|
||||
call ijk2idx(idx,[i,j,k],[nx,ny,nz],base)
|
||||
end subroutine lijk2lidx3d
|
||||
|
||||
subroutine lijk2lidx2d(idx,i,j,nx,ny,base)
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_epk_), intent(out) :: idx
|
||||
integer(psb_epk_), intent(in) :: i,j,nx,ny
|
||||
integer(psb_mpk_), intent(in), optional :: base
|
||||
|
||||
! idx = ((i-base_)*ny + (j-base_) + base_
|
||||
call ijk2idx(idx,[i,j],[nx,ny],base)
|
||||
end subroutine lijk2lidx2d
|
||||
|
||||
!
|
||||
! dist1Didx
|
||||
! Given an index space [base : N-(1-base)] and
|
||||
@@ -263,7 +592,7 @@ contains
|
||||
! V(P+1) : first index owned by process P+1
|
||||
!
|
||||
subroutine dist1Didx(v,n,np,base)
|
||||
use psb_base_mod, only : psb_ipk_
|
||||
use psb_base_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(out) :: v(:)
|
||||
integer(psb_ipk_), intent(in) :: n, np
|
||||
|
||||
@@ -314,3 +314,290 @@ subroutine shb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine shb_write
|
||||
|
||||
subroutine lshb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
real(psb_spk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
|
||||
character :: rhstype*3,type*3,key*8
|
||||
character(len=72) :: mtitle_
|
||||
character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
|
||||
integer(psb_lpk_) :: indcrd, ptrcrd, totcrd,&
|
||||
& valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
|
||||
type(psb_ls_csc_sparse_mat) :: acsc
|
||||
type(psb_ls_coo_sparse_mat) :: acoo
|
||||
integer(psb_ipk_) :: ircode, infile, info
|
||||
integer(psb_lpk_) :: i,nzr
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
iret = 0
|
||||
ircode = 0
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
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 acsc%allocate(nrow,ncol,nnzero)
|
||||
if (ircode /= 0 ) then
|
||||
write(psb_err_unit,*) 'Memory allocation failed'
|
||||
goto 993
|
||||
end if
|
||||
|
||||
if (present(mtitle)) mtitle=mtitle_
|
||||
|
||||
|
||||
if (psb_tolower(type(1:1)) == 'r') then
|
||||
if (psb_tolower(type(2:2)) == 'u') then
|
||||
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
call a%mv_from(acsc)
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 's') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
|
||||
call a%cscnv(ircode,type='csr')
|
||||
if (infile /= 5) close(infile)
|
||||
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 iret=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 iret=902
|
||||
write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
|
||||
return
|
||||
993 iret=993
|
||||
write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
|
||||
return
|
||||
end subroutine lshb_read
|
||||
|
||||
subroutine lshb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(in), target :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
real(psb_spk_), optional :: rhs(:), g(:), x(:)
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
character(len=*), parameter:: ptrfmt='(10I8)',indfmt='(10I8)'
|
||||
integer(psb_ipk_), parameter :: jptr=10,jind=10
|
||||
character(len=*), parameter:: valfmt='(4E20.12)',rhsfmt='(4E20.12)'
|
||||
integer(psb_ipk_), parameter :: jval=4,jrhs=4
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
type(psb_ls_csc_sparse_mat) :: acsc
|
||||
character(len=72) :: mtitle_
|
||||
character(len=8) :: key_
|
||||
|
||||
character :: rhstype*3,type*3
|
||||
|
||||
integer(psb_ipk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
|
||||
& nrow,ncol,nnzero, neltvl, nrhs, nrhsix
|
||||
|
||||
iret = 0
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
if (present(mtitle)) then
|
||||
mtitle_ = mtitle
|
||||
else
|
||||
mtitle_ = 'Temporary PSBLAS title '
|
||||
endif
|
||||
if (present(key)) then
|
||||
key_ = key
|
||||
else
|
||||
key_ = 'PSBMAT00'
|
||||
endif
|
||||
|
||||
|
||||
call acsc%cp_from_fmt(a%a, iret)
|
||||
if (iret /= 0) return
|
||||
|
||||
|
||||
nrow = acsc%get_nrows()
|
||||
ncol = acsc%get_ncols()
|
||||
nnzero = acsc%get_nzeros()
|
||||
|
||||
neltvl = 0
|
||||
|
||||
ptrcrd = (ncol+1)/jptr
|
||||
if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
|
||||
indcrd = nnzero/jind
|
||||
if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
|
||||
valcrd = nnzero/jval
|
||||
if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
|
||||
rhstype = ''
|
||||
if (present(rhs)) then
|
||||
if (size(rhs)<nrow) then
|
||||
rhscrd = 0
|
||||
else
|
||||
rhscrd = nrow/jrhs
|
||||
if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
|
||||
endif
|
||||
nrhs = 1
|
||||
rhstype(1:1) = 'F'
|
||||
else
|
||||
rhscrd = 0
|
||||
nrhs = 0
|
||||
end if
|
||||
totcrd = ptrcrd + indcrd + valcrd + rhscrd
|
||||
|
||||
nrhsix = nrhs*nrow
|
||||
|
||||
if (present(g)) then
|
||||
rhstype(2:2) = 'G'
|
||||
end if
|
||||
if (present(x)) then
|
||||
rhstype(3:3) = 'X'
|
||||
end if
|
||||
type = 'RUA'
|
||||
|
||||
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) (acsc%icp(i),i=1,ncol+1)
|
||||
write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) write (iout,fmt=valfmt) (acsc%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)
|
||||
|
||||
|
||||
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
iret=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine lshb_write
|
||||
|
||||
@@ -88,11 +88,11 @@ subroutine psb_smatdist(a_glob, a, ictxt, desc_a,&
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing
|
||||
integer(psb_ipk_) :: i_count, j_count,&
|
||||
& k_count, root, liwork, nrow, ncol, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act, int_err(5)
|
||||
integer(psb_ipk_) :: k_count, root, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_ipk_), allocatable :: irow(:),icol(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
real(psb_spk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
@@ -140,8 +140,7 @@ subroutine psb_smatdist(a_glob, a, ictxt, desc_a,&
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
int_err(1)=liwork
|
||||
call psb_errpush(info,name,i_err=int_err,a_err='integer')
|
||||
call psb_errpush(info,name,i_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
@@ -339,3 +338,319 @@ subroutine psb_smatdist(a_glob, a, ictxt, desc_a,&
|
||||
return
|
||||
|
||||
end subroutine psb_smatdist
|
||||
|
||||
|
||||
|
||||
subroutine psb_lsmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_sspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_sspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
use psb_base_mod
|
||||
use psb_mat_mod
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lsspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_sspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_s_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing, root, iproc
|
||||
integer(psb_ipk_) :: err_act, il, inz
|
||||
integer(psb_lpk_) :: k_count, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, nnr, err
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
real(psb_spk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info = psb_success_
|
||||
err = 0
|
||||
name = 'psb_s_mat_dist'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
! executable statements
|
||||
if (present(inroot)) then
|
||||
root = inroot
|
||||
else
|
||||
root = psb_root_
|
||||
end if
|
||||
call psb_info(ictxt, iam, np)
|
||||
if (iam == root) then
|
||||
nrow = a_glob%get_nrows()
|
||||
ncol = a_glob%get_ncols()
|
||||
if (nrow /= ncol) then
|
||||
write(psb_err_unit,*) 'a rectangular matrix ? ',nrow,ncol
|
||||
info=-1
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
endif
|
||||
nnzero = a_glob%get_nzeros()
|
||||
nrhs = 1
|
||||
endif
|
||||
|
||||
use_parts = present(parts)
|
||||
use_v = present(v)
|
||||
if (count((/ use_parts, use_v /)) /= 1) then
|
||||
info=psb_err_no_optional_arg_
|
||||
call psb_errpush(info,name,a_err=" v, parts")
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
! broadcast informations to other processors
|
||||
call psb_bcast(ictxt,nrow, root)
|
||||
call psb_bcast(ictxt,ncol, root)
|
||||
call psb_bcast(ictxt,nnzero, root)
|
||||
call psb_bcast(ictxt,nrhs, root)
|
||||
liwork = max(np, nrow + ncol)
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,l_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
write (*, fmt = *) 'start matdist',root, size(iwork),&
|
||||
&nrow, ncol, nnzero,nrhs
|
||||
endif
|
||||
if (use_parts) then
|
||||
call psb_cdall(ictxt,desc_a,info,mg=nrow,parts=parts)
|
||||
else if (use_v) then
|
||||
call psb_cdall(ictxt,desc_a,info,vg=v)
|
||||
else
|
||||
info = -1
|
||||
end if
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
inz = ((nnzero+np-1)/np)
|
||||
call psb_spall(a,desc_a,info,nnz=inz)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
isize = 3*nb*max(((nnzero+nrow)/nrow),nb)
|
||||
allocate(val(isize),irow(isize),icol(isize),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
i_count = 1
|
||||
|
||||
do while (i_count <= nrow)
|
||||
|
||||
if (use_parts) then
|
||||
call parts(i_count,nrow,np,iwork, np_sharing)
|
||||
!
|
||||
! np_sharing allows for overlap in the data distribution.
|
||||
! If an index is overlapped, then we have to send its row
|
||||
! to multiple processes. NOTE: we are assuming the output
|
||||
! from PARTS is coherent, otherwise a deadlock is the most
|
||||
! likely outcome.
|
||||
!
|
||||
j_count = i_count
|
||||
if (np_sharing == 1) then
|
||||
iproc = iwork(1)
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
call parts(j_count,nrow,np,iwrk2, np_sharing)
|
||||
if (np_sharing /= 1 ) exit
|
||||
if (iwrk2(1) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
else
|
||||
np_sharing = 1
|
||||
j_count = i_count
|
||||
iproc = v(i_count)
|
||||
iwork(1:np_sharing) = iproc
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
if (v(j_count) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
|
||||
! now we should insert rows i_count..j_count-1
|
||||
nnr = j_count - i_count
|
||||
|
||||
if (iam == root) then
|
||||
|
||||
ll = 0
|
||||
do i= i_count, j_count-1
|
||||
call a_glob%csget(i,i,nz,&
|
||||
& irow,icol,val,info,nzin=ll,append=.true.)
|
||||
if (info /= psb_success_) then
|
||||
if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then
|
||||
write(psb_err_unit,*) 'Allocation failure? This should not happen!'
|
||||
end if
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
ll = ll + nz
|
||||
end do
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
|
||||
if (iproc == iam) then
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
else
|
||||
call psb_snd(ictxt,nnr,iproc)
|
||||
call psb_snd(ictxt,ll,iproc)
|
||||
call psb_snd(ictxt,irow(1:ll),iproc)
|
||||
call psb_snd(ictxt,icol(1:ll),iproc)
|
||||
call psb_snd(ictxt,val(1:ll),iproc)
|
||||
call psb_rcv(ictxt,ll,iproc)
|
||||
endif
|
||||
end do
|
||||
else if (iam /= root) then
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
if (iproc == iam) then
|
||||
call psb_rcv(ictxt,nnr,root)
|
||||
call psb_rcv(ictxt,ll,root)
|
||||
if (ll > size(irow)) then
|
||||
write(psb_err_unit,*) iam,'need to reallocate ',ll
|
||||
deallocate(val,irow,icol)
|
||||
allocate(val(ll),irow(ll),icol(ll),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
endif
|
||||
call psb_rcv(ictxt,irow(1:ll),root)
|
||||
call psb_rcv(ictxt,icol(1:ll),root)
|
||||
call psb_rcv(ictxt,val(1:ll),root)
|
||||
call psb_snd(ictxt,ll,root)
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psspins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
end do
|
||||
endif
|
||||
i_count = j_count
|
||||
end do
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t0 = psb_wtime()
|
||||
call psb_cdasb(desc_a,info)
|
||||
t1 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime()
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=fmt,mold=mold)
|
||||
t3 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (iam == root) then
|
||||
write(psb_out_unit,*) 'descriptor assembly: ',t1-t0
|
||||
write(psb_out_unit,*) 'sparse matrix assembly: ',t3-t2
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(val,irow,icol,iwork,iwrk2,stat=info)
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='deallocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (iam == root) write (*, fmt = *) 'end matdist'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ictxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
end subroutine psb_lsmatdist
|
||||
|
||||
@@ -31,7 +31,8 @@
|
||||
!
|
||||
module psb_s_mat_dist_mod
|
||||
use psb_base_mod, only : psb_ipk_, psb_spk_, psb_desc_type, psb_parts, &
|
||||
& psb_sspmat_type, psb_s_base_sparse_mat, psb_s_vect_type
|
||||
& psb_sspmat_type, psb_s_base_sparse_mat, psb_s_vect_type, &
|
||||
& psb_lsspmat_type
|
||||
|
||||
interface psb_matdist
|
||||
subroutine psb_smatdist(a_glob, a, ictxt, desc_a,&
|
||||
@@ -91,6 +92,64 @@ module psb_s_mat_dist_mod
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_smatdist
|
||||
subroutine psb_lsmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_lsspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_sspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
import :: psb_ipk_, psb_sspmat_type, psb_spk_, psb_desc_type,&
|
||||
& psb_s_base_sparse_mat, psb_s_vect_type, psb_parts, &
|
||||
& psb_lsspmat_type
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lsspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_sspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_s_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_lsmatdist
|
||||
end interface
|
||||
|
||||
end module psb_s_mat_dist_mod
|
||||
|
||||
@@ -434,3 +434,156 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename)
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine smm_mat_write
|
||||
|
||||
subroutine lsmm_mat_read(a, info, iunit, filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character :: mmheader*15, fmt*15, object*10, type*10, sym*15
|
||||
character(1024) :: line
|
||||
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr
|
||||
integer(psb_ipk_) :: ircode,infile
|
||||
type(psb_ls_coo_sparse_mat), allocatable :: acoo
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
read(infile,fmt=*,end=902) mmheader, object, fmt, type, sym
|
||||
|
||||
if ( (psb_tolower(object) /= 'matrix').or.(psb_tolower(fmt) /= 'coordinate')) then
|
||||
write(psb_err_unit,*) 'READ_MATRIX: input file type not yet supported'
|
||||
info=909
|
||||
return
|
||||
end if
|
||||
|
||||
do
|
||||
read(infile,fmt='(a)') line
|
||||
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 acoo%allocate(nrow,ncol,nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
|
||||
end do
|
||||
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 acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(info)
|
||||
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
info=904
|
||||
end if
|
||||
|
||||
|
||||
if (infile /= 5) close(infile)
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 info=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 info=902
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Unexpected end of file '
|
||||
return
|
||||
905 info=905
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Error at line',i
|
||||
return
|
||||
993 info=993
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Memory allocation failure'
|
||||
return
|
||||
end subroutine lsmm_mat_read
|
||||
|
||||
|
||||
subroutine lsmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lsspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
call a%print(iout,head=mtitle)
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
info=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine lsmm_mat_write
|
||||
|
||||
@@ -268,7 +268,7 @@ contains
|
||||
name = 'mat_renum_amd'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
#if defined(HAVE_AMD)
|
||||
#if defined(HAVE_AMD) && defined(IPK4)
|
||||
|
||||
info = psb_success_
|
||||
nr = a%get_nrows()
|
||||
|
||||
+331
-1
@@ -349,6 +349,336 @@ subroutine zhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
|
||||
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
iret=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine zhb_write
|
||||
|
||||
subroutine lzhb_read(a, iret, iunit, filename,b,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
complex(psb_dpk_), optional, allocatable, intent(out) :: b(:,:), g(:,:), x(:,:)
|
||||
character(len=72), optional, intent(out) :: mtitle
|
||||
|
||||
character :: rhstype*3,type*3,key*8
|
||||
character(len=72) :: mtitle_
|
||||
character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
|
||||
integer(psb_lpk_) :: indcrd, ptrcrd, totcrd,&
|
||||
& valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
|
||||
type(psb_lz_csc_sparse_mat) :: acsc
|
||||
type(psb_lz_coo_sparse_mat) :: acoo
|
||||
integer(psb_ipk_) :: ircode, infile, info
|
||||
integer(psb_lpk_) :: i,nzr
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
iret = 0
|
||||
ircode = 0
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
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 acsc%allocate(nrow,ncol,nnzero)
|
||||
if (ircode /= 0 ) then
|
||||
write(psb_err_unit,*) 'Memory allocation failed'
|
||||
goto 993
|
||||
end if
|
||||
|
||||
if (present(mtitle)) mtitle=mtitle_
|
||||
|
||||
|
||||
if (psb_tolower(type(1:1)) == 'c') then
|
||||
if (psb_tolower(type(2:2)) == 'u') then
|
||||
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
call a%mv_from(acsc)
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 's') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else if (psb_tolower(type(2:2)) == 'h') then
|
||||
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
|
||||
read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
|
||||
read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
|
||||
|
||||
|
||||
if (present(b)) then
|
||||
if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,b,info)
|
||||
read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(g)) then
|
||||
if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,g,info)
|
||||
read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
if (present(x)) then
|
||||
if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then
|
||||
call psb_realloc(nrow,1,x,info)
|
||||
read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
call acoo%mv_from_fmt(acsc,info)
|
||||
call acoo%reallocate(2*nnzero)
|
||||
! A is now in COO format
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = conjg(acoo%val(i))
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(ircode)
|
||||
if (ircode == 0) call a%mv_from(acoo)
|
||||
if (ircode /= 0) goto 993
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
iret=904
|
||||
end if
|
||||
|
||||
call a%cscnv(ircode,type='csr')
|
||||
if (infile /= 5) close(infile)
|
||||
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 iret=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 iret=902
|
||||
write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
|
||||
return
|
||||
993 iret=993
|
||||
write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
|
||||
return
|
||||
end subroutine lzhb_read
|
||||
|
||||
subroutine lzhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(in), target :: a
|
||||
integer(psb_ipk_), intent(out) :: iret
|
||||
character(len=*), optional, intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character(len=*), optional, intent(in) :: key
|
||||
complex(psb_dpk_), optional :: rhs(:), g(:), x(:)
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
character(len=*), parameter:: ptrfmt='(10I8)',indfmt='(10I8)'
|
||||
integer(psb_ipk_), parameter :: jptr=10,jind=10
|
||||
character(len=*), parameter:: valfmt='(4E20.12)',rhsfmt='(4E20.12)'
|
||||
integer(psb_ipk_), parameter :: jval=2,jrhs=2
|
||||
character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
|
||||
character(len=*), parameter :: fmt11='(a3,11x,2i14)'
|
||||
character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
|
||||
|
||||
type(psb_lz_csc_sparse_mat) :: acsc
|
||||
character(len=72) :: mtitle_
|
||||
character(len=8) :: key_
|
||||
|
||||
character :: rhstype*3,type*3
|
||||
|
||||
integer(psb_ipk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
|
||||
& nrow,ncol,nnzero, neltvl, nrhs, nrhsix
|
||||
|
||||
iret = 0
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
if (present(mtitle)) then
|
||||
mtitle_ = mtitle
|
||||
else
|
||||
mtitle_ = 'Temporary PSBLAS title '
|
||||
endif
|
||||
if (present(key)) then
|
||||
key_ = key
|
||||
else
|
||||
key_ = 'PSBMAT00'
|
||||
endif
|
||||
|
||||
|
||||
call acsc%cp_from_fmt(a%a, iret)
|
||||
if (iret /= 0) return
|
||||
|
||||
|
||||
nrow = acsc%get_nrows()
|
||||
ncol = acsc%get_ncols()
|
||||
nnzero = acsc%get_nzeros()
|
||||
|
||||
neltvl = 0
|
||||
|
||||
ptrcrd = (ncol+1)/jptr
|
||||
if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
|
||||
indcrd = nnzero/jind
|
||||
if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
|
||||
valcrd = nnzero/jval
|
||||
if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
|
||||
rhstype = ''
|
||||
if (present(rhs)) then
|
||||
if (size(rhs)<nrow) then
|
||||
rhscrd = 0
|
||||
else
|
||||
rhscrd = nrow/jrhs
|
||||
if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
|
||||
endif
|
||||
nrhs = 1
|
||||
rhstype(1:1) = 'F'
|
||||
else
|
||||
rhscrd = 0
|
||||
nrhs = 0
|
||||
end if
|
||||
totcrd = ptrcrd + indcrd + valcrd + rhscrd
|
||||
|
||||
nrhsix = nrhs*nrow
|
||||
|
||||
if (present(g)) then
|
||||
rhstype(2:2) = 'G'
|
||||
end if
|
||||
if (present(x)) then
|
||||
rhstype(3:3) = 'X'
|
||||
end if
|
||||
type = 'CUA'
|
||||
|
||||
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) (acsc%icp(i),i=1,ncol+1)
|
||||
write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
|
||||
if (valcrd > 0) write (iout,fmt=valfmt) (acsc%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)
|
||||
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
@@ -360,4 +690,4 @@ subroutine zhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
|
||||
iret=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine zhb_write
|
||||
end subroutine lzhb_write
|
||||
|
||||
@@ -88,11 +88,11 @@ subroutine psb_zmatdist(a_glob, a, ictxt, desc_a,&
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing
|
||||
integer(psb_ipk_) :: i_count, j_count,&
|
||||
& k_count, root, liwork, nrow, ncol, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act, int_err(5)
|
||||
integer(psb_ipk_) :: k_count, root, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, iproc, nnr, err, err_act
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_ipk_), allocatable :: irow(:),icol(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
complex(psb_dpk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
@@ -140,8 +140,7 @@ subroutine psb_zmatdist(a_glob, a, ictxt, desc_a,&
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
int_err(1)=liwork
|
||||
call psb_errpush(info,name,i_err=int_err,a_err='integer')
|
||||
call psb_errpush(info,name,i_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
@@ -339,3 +338,319 @@ subroutine psb_zmatdist(a_glob, a, ictxt, desc_a,&
|
||||
return
|
||||
|
||||
end subroutine psb_zmatdist
|
||||
|
||||
|
||||
|
||||
subroutine psb_lzmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_zspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_zspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
use psb_base_mod
|
||||
use psb_mat_mod
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lzspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_zspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_z_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
|
||||
! local variables
|
||||
logical :: use_parts, use_v
|
||||
integer(psb_ipk_) :: np, iam, np_sharing, root, iproc
|
||||
integer(psb_ipk_) :: err_act, il, inz
|
||||
integer(psb_lpk_) :: k_count, liwork, nnzero, nrhs,&
|
||||
& i, ll, nz, isize, nnr, err
|
||||
integer(psb_lpk_) :: i_count, j_count, nrow, ncol, ig
|
||||
integer(psb_ipk_), allocatable :: iwork(:), iwrk2(:)
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:)
|
||||
complex(psb_dpk_), allocatable :: val(:)
|
||||
integer(psb_ipk_), parameter :: nb=30
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, t4, t5
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info = psb_success_
|
||||
err = 0
|
||||
name = 'psb_z_mat_dist'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
! executable statements
|
||||
if (present(inroot)) then
|
||||
root = inroot
|
||||
else
|
||||
root = psb_root_
|
||||
end if
|
||||
call psb_info(ictxt, iam, np)
|
||||
if (iam == root) then
|
||||
nrow = a_glob%get_nrows()
|
||||
ncol = a_glob%get_ncols()
|
||||
if (nrow /= ncol) then
|
||||
write(psb_err_unit,*) 'a rectangular matrix ? ',nrow,ncol
|
||||
info=-1
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
endif
|
||||
nnzero = a_glob%get_nzeros()
|
||||
nrhs = 1
|
||||
endif
|
||||
|
||||
use_parts = present(parts)
|
||||
use_v = present(v)
|
||||
if (count((/ use_parts, use_v /)) /= 1) then
|
||||
info=psb_err_no_optional_arg_
|
||||
call psb_errpush(info,name,a_err=" v, parts")
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
! broadcast informations to other processors
|
||||
call psb_bcast(ictxt,nrow, root)
|
||||
call psb_bcast(ictxt,ncol, root)
|
||||
call psb_bcast(ictxt,nnzero, root)
|
||||
call psb_bcast(ictxt,nrhs, root)
|
||||
liwork = max(np, nrow + ncol)
|
||||
allocate(iwork(liwork), iwrk2(np),stat = info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,l_err=(/liwork/),a_err='integer')
|
||||
goto 9999
|
||||
endif
|
||||
if (iam == root) then
|
||||
write (*, fmt = *) 'start matdist',root, size(iwork),&
|
||||
&nrow, ncol, nnzero,nrhs
|
||||
endif
|
||||
if (use_parts) then
|
||||
call psb_cdall(ictxt,desc_a,info,mg=nrow,parts=parts)
|
||||
else if (use_v) then
|
||||
call psb_cdall(ictxt,desc_a,info,vg=v)
|
||||
else
|
||||
info = -1
|
||||
end if
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
inz = ((nnzero+np-1)/np)
|
||||
call psb_spall(a,desc_a,info,nnz=inz)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spall'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
isize = 3*nb*max(((nnzero+nrow)/nrow),nb)
|
||||
allocate(val(isize),irow(isize),icol(isize),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
i_count = 1
|
||||
|
||||
do while (i_count <= nrow)
|
||||
|
||||
if (use_parts) then
|
||||
call parts(i_count,nrow,np,iwork, np_sharing)
|
||||
!
|
||||
! np_sharing allows for overlap in the data distribution.
|
||||
! If an index is overlapped, then we have to send its row
|
||||
! to multiple processes. NOTE: we are assuming the output
|
||||
! from PARTS is coherent, otherwise a deadlock is the most
|
||||
! likely outcome.
|
||||
!
|
||||
j_count = i_count
|
||||
if (np_sharing == 1) then
|
||||
iproc = iwork(1)
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
call parts(j_count,nrow,np,iwrk2, np_sharing)
|
||||
if (np_sharing /= 1 ) exit
|
||||
if (iwrk2(1) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
else
|
||||
np_sharing = 1
|
||||
j_count = i_count
|
||||
iproc = v(i_count)
|
||||
iwork(1:np_sharing) = iproc
|
||||
do
|
||||
j_count = j_count + 1
|
||||
if (j_count-i_count >= nb) exit
|
||||
if (j_count > nrow) exit
|
||||
if (v(j_count) /= iproc ) exit
|
||||
end do
|
||||
end if
|
||||
|
||||
! now we should insert rows i_count..j_count-1
|
||||
nnr = j_count - i_count
|
||||
|
||||
if (iam == root) then
|
||||
|
||||
ll = 0
|
||||
do i= i_count, j_count-1
|
||||
call a_glob%csget(i,i,nz,&
|
||||
& irow,icol,val,info,nzin=ll,append=.true.)
|
||||
if (info /= psb_success_) then
|
||||
if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then
|
||||
write(psb_err_unit,*) 'Allocation failure? This should not happen!'
|
||||
end if
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
ll = ll + nz
|
||||
end do
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
|
||||
if (iproc == iam) then
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
else
|
||||
call psb_snd(ictxt,nnr,iproc)
|
||||
call psb_snd(ictxt,ll,iproc)
|
||||
call psb_snd(ictxt,irow(1:ll),iproc)
|
||||
call psb_snd(ictxt,icol(1:ll),iproc)
|
||||
call psb_snd(ictxt,val(1:ll),iproc)
|
||||
call psb_rcv(ictxt,ll,iproc)
|
||||
endif
|
||||
end do
|
||||
else if (iam /= root) then
|
||||
|
||||
do k_count = 1, np_sharing
|
||||
iproc = iwork(k_count)
|
||||
if (iproc == iam) then
|
||||
call psb_rcv(ictxt,nnr,root)
|
||||
call psb_rcv(ictxt,ll,root)
|
||||
if (ll > size(irow)) then
|
||||
write(psb_err_unit,*) iam,'need to reallocate ',ll
|
||||
deallocate(val,irow,icol)
|
||||
allocate(val(ll),irow(ll),icol(ll),stat=info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='Allocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
endif
|
||||
call psb_rcv(ictxt,irow(1:ll),root)
|
||||
call psb_rcv(ictxt,icol(1:ll),root)
|
||||
call psb_rcv(ictxt,val(1:ll),root)
|
||||
call psb_snd(ictxt,ll,root)
|
||||
il = ll
|
||||
call psb_spins(il,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psspins'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
end do
|
||||
endif
|
||||
i_count = j_count
|
||||
end do
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t0 = psb_wtime()
|
||||
call psb_cdasb(desc_a,info)
|
||||
t1 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_cdasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_barrier(ictxt)
|
||||
t2 = psb_wtime()
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=fmt,mold=mold)
|
||||
t3 = psb_wtime()
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_spasb'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
if (iam == root) then
|
||||
write(psb_out_unit,*) 'descriptor assembly: ',t1-t0
|
||||
write(psb_out_unit,*) 'sparse matrix assembly: ',t3-t2
|
||||
end if
|
||||
|
||||
|
||||
|
||||
deallocate(val,irow,icol,iwork,iwrk2,stat=info)
|
||||
if(info /= psb_success_)then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='deallocate'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
if (iam == root) write (*, fmt = *) 'end matdist'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ictxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
end subroutine psb_lzmatdist
|
||||
|
||||
@@ -31,7 +31,8 @@
|
||||
!
|
||||
module psb_z_mat_dist_mod
|
||||
use psb_base_mod, only : psb_ipk_, psb_dpk_, psb_desc_type, psb_parts, &
|
||||
& psb_zspmat_type, psb_z_base_sparse_mat, psb_z_vect_type
|
||||
& psb_zspmat_type, psb_z_base_sparse_mat, psb_z_vect_type, &
|
||||
& psb_lzspmat_type
|
||||
|
||||
interface psb_matdist
|
||||
subroutine psb_zmatdist(a_glob, a, ictxt, desc_a,&
|
||||
@@ -91,6 +92,64 @@ module psb_z_mat_dist_mod
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_zmatdist
|
||||
subroutine psb_lzmatdist(a_glob, a, ictxt, desc_a,&
|
||||
& info, parts, v, inroot,fmt,mold)
|
||||
!
|
||||
! an utility subroutine to distribute a matrix among processors
|
||||
! according to a user defined data distribution, using
|
||||
! sparse matrix subroutines.
|
||||
!
|
||||
! type(psb_lzspmat) :: a_glob
|
||||
! on entry: this contains the global sparse matrix as follows:
|
||||
!
|
||||
! type(psb_zspmat_type) :: a
|
||||
! on exit : this will contain the local sparse matrix.
|
||||
!
|
||||
! interface parts
|
||||
! ! .....user passed subroutine.....
|
||||
! subroutine parts(global_indx,n,np,pv,nv)
|
||||
! implicit none
|
||||
! integer(psb_ipk_), intent(in) :: global_indx, n, np
|
||||
! integer(psb_ipk_), intent(out) :: nv
|
||||
! integer(psb_ipk_), intent(out) :: pv(*)
|
||||
!
|
||||
! end subroutine parts
|
||||
! end interface
|
||||
! on entry: subroutine providing user defined data distribution.
|
||||
! for each global_indx the subroutine should return
|
||||
! the list pv of all processes owning the row with
|
||||
! that index; the list will contain nv entries.
|
||||
! usually nv=1; if nv >1 then we have an overlap in the data
|
||||
! distribution.
|
||||
!
|
||||
! integer(psb_ipk_) :: ictxt
|
||||
! on entry: the PSBLAS parallel environment context.
|
||||
!
|
||||
! type (desc_type) :: desc_a
|
||||
! on exit : the updated array descriptor
|
||||
!
|
||||
!
|
||||
! integer(psb_ipk_), optional :: inroot
|
||||
! on entry: specifies processor holding a_glob. default: 0
|
||||
! on exit : unchanged.
|
||||
!
|
||||
import :: psb_ipk_, psb_zspmat_type, psb_dpk_, psb_desc_type,&
|
||||
& psb_z_base_sparse_mat, psb_z_vect_type, psb_parts, &
|
||||
& psb_lzspmat_type
|
||||
implicit none
|
||||
|
||||
! parameters
|
||||
type(psb_lzspmat_type) :: a_glob
|
||||
integer(psb_ipk_) :: ictxt
|
||||
type(psb_zspmat_type) :: a
|
||||
type(psb_desc_type) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional :: inroot
|
||||
character(len=*), optional :: fmt
|
||||
class(psb_z_base_sparse_mat), optional :: mold
|
||||
procedure(psb_parts), optional :: parts
|
||||
integer(psb_ipk_), optional :: v(:)
|
||||
end subroutine psb_lzmatdist
|
||||
end interface
|
||||
|
||||
end module psb_z_mat_dist_mod
|
||||
|
||||
@@ -462,3 +462,178 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
end subroutine zmm_mat_write
|
||||
|
||||
|
||||
subroutine lzmm_mat_read(a, info, iunit, filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(out) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
character :: mmheader*15, fmt*15, object*10, type*10, sym*15
|
||||
character(1024) :: line
|
||||
integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr
|
||||
integer(psb_ipk_) :: ircode,infile
|
||||
type(psb_lz_coo_sparse_mat), allocatable :: acoo
|
||||
real(psb_dpk_) :: are, aim
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
infile=5
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=99
|
||||
endif
|
||||
open(infile,file=filename, status='OLD', err=901, action='READ')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
infile=iunit
|
||||
else
|
||||
infile=5
|
||||
endif
|
||||
endif
|
||||
|
||||
read(infile,fmt=*,end=902) mmheader, object, fmt, type, sym
|
||||
|
||||
if ( (psb_tolower(object) /= 'matrix').or.(psb_tolower(fmt) /= 'coordinate')) then
|
||||
write(psb_err_unit,*) 'READ_MATRIX: input file type not yet supported'
|
||||
info=909
|
||||
return
|
||||
end if
|
||||
|
||||
do
|
||||
read(infile,fmt='(a)') line
|
||||
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) == 'complex').and.(psb_tolower(sym) == 'general')) then
|
||||
call acoo%allocate(nrow,ncol,nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
|
||||
end do
|
||||
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) == 'complex').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 acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = acoo%val(i)
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(info)
|
||||
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
|
||||
else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then
|
||||
! we are generally working with non-symmetric matrices, so
|
||||
! we de-symmetrize what we are about to read
|
||||
call acoo%allocate(nrow,ncol,2*nnzero)
|
||||
do i=1,nnzero
|
||||
read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim
|
||||
acoo%val(i) = cmplx(are,aim,kind=psb_dpk_)
|
||||
end do
|
||||
nzr = nnzero
|
||||
do i=1,nnzero
|
||||
if (acoo%ia(i) /= acoo%ja(i)) then
|
||||
nzr = nzr + 1
|
||||
acoo%val(nzr) = conjg(acoo%val(i))
|
||||
acoo%ia(nzr) = acoo%ja(i)
|
||||
acoo%ja(nzr) = acoo%ia(i)
|
||||
end if
|
||||
end do
|
||||
call acoo%set_nzeros(nzr)
|
||||
call acoo%fix(info)
|
||||
|
||||
call a%mv_from(acoo)
|
||||
call a%cscnv(ircode,type='csr')
|
||||
|
||||
else
|
||||
write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
|
||||
info=904
|
||||
end if
|
||||
if (infile /= 5) close(infile)
|
||||
return
|
||||
|
||||
! open failed
|
||||
901 info=901
|
||||
write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
|
||||
return
|
||||
902 info=902
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Unexpected end of file '
|
||||
return
|
||||
905 info=905
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Error at line',i
|
||||
return
|
||||
993 info=993
|
||||
write(psb_err_unit,*) 'READ_MATRIX: Memory allocation failure'
|
||||
return
|
||||
end subroutine lzmm_mat_read
|
||||
|
||||
subroutine lzmm_mat_write(a,mtitle,info,iunit,filename)
|
||||
use psb_base_mod
|
||||
implicit none
|
||||
type(psb_lzspmat_type), intent(in) :: a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character(len=*), intent(in) :: mtitle
|
||||
integer(psb_ipk_), optional, intent(in) :: iunit
|
||||
character(len=*), optional, intent(in) :: filename
|
||||
integer(psb_ipk_) :: iout
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
if (present(filename)) then
|
||||
if (filename == '-') then
|
||||
iout=6
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=99
|
||||
endif
|
||||
open(iout,file=filename, err=901, action='WRITE')
|
||||
endif
|
||||
else
|
||||
if (present(iunit)) then
|
||||
iout = iunit
|
||||
else
|
||||
iout=6
|
||||
endif
|
||||
endif
|
||||
|
||||
call a%print(iout,head=mtitle)
|
||||
|
||||
if (iout /= 6) close(iout)
|
||||
|
||||
|
||||
return
|
||||
|
||||
901 continue
|
||||
info=901
|
||||
write(psb_err_unit,*) 'Error while opening ',filename
|
||||
return
|
||||
end subroutine lzmm_mat_write
|
||||
|
||||
|
||||
@@ -267,7 +267,7 @@ contains
|
||||
name = 'mat_renum_amd'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
#if defined(HAVE_AMD)
|
||||
#if defined(HAVE_AMD) && defined(IPK4)
|
||||
|
||||
info = psb_success_
|
||||
nr = a%get_nrows()
|
||||
|
||||
Reference in New Issue
Block a user