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

psblas3:

Browse Source 
base/comm/Makefile
 base/comm/psb_dspgather.F90
 base/modules/Makefile
 base/modules/psb_comm_mod.f90
 base/modules/psb_d_mat_mod.f03
 base/modules/psb_serial_mod.f90
 base/serial/Makefile
 base/serial/f77/Makefile
 base/serial/psb_dnumbmm.f90
 base/serial/psb_drwextd.f90
 base/serial/psb_dsymbmm.f90
 configure.ac
 configure

Beginning work on new MLD, need support for other operators.
This commit is contained in:
Salvatore Filippone
2009-10-09 08:36:59 +00:00
parent ebe0d004ea
commit eb0896c378
13 changed files with 838 additions and 263 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/comm/Makefile
+3 -1
View File
@@ -6,7 +6,9 @@ OBJS = psb_dgather.o psb_dhalo.o psb_dovrl.o \
psb_cgather.o psb_chalo.o psb_covrl.o \
psb_zgather.o psb_zhalo.o psb_zovrl.o
MPFOBJS=psb_dscatter.o psb_zscatter.o psb_iscatter.o psb_cscatter.o psb_sscatter.o
MPFOBJS=psb_dscatter.o psb_zscatter.o psb_iscatter.o psb_cscatter.o psb_sscatter.o\
psb_dspgather.o
LIBDIR=..
MODDIR=../modules
FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG)$(MODDIR) $(FMFLAG).
base/comm/psb_dspgather.F90
+113
View File
@@ -0,0 +1,113 @@
subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
use psb_descriptor_type
use psb_error_mod
use psb_mat_mod
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
type(psb_d_sparse_mat), intent(inout) :: loca
type(psb_d_sparse_mat), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer, intent(in), optional :: root, dupl
logical, intent(in), optional :: keepnum,keeploc
type(psb_d_coo_sparse_mat) :: loc_coo, glob_coo
integer :: ictxt,np,me, err_act, icomm, dupl_, nrg, ncg, nzg
integer :: ip, ndx,naggrm1,naggrp1, i, j, k
logical :: keepnum_, keeploc_
integer, allocatable :: nzbr(:), idisp(:)
character(len=20) :: name
integer :: debug_level, debug_unit
name='psb_gather'
if (psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_a)
icomm = psb_cd_get_mpic(desc_a)
call psb_info(ictxt, me, np)
if (present(keepnum)) then
keepnum_ = keepnum
else
keepnum_ = .true.
end if
if (present(keeploc)) then
keeploc_ = keeploc
else
keeploc_ = .true.
end if
if (keepnum_) then
nrg = psb_cd_get_global_rows(desc_a)
ncg = psb_cd_get_global_rows(desc_a)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
goto 9999
end if
call loca%mv_to(loc_coo)
nzbr(:) = 0
nzbr(me+1) = loc_coo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzg = sum(nzbr)
if (info == 0) call glob_coo%allocate(nrg,ncg,nzg)
if (info /= 0) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(loc_coo%val,ndx,mpi_double_precision,&
& glob_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == 0) call mpi_allgatherv(loc_coo%ia,ndx,mpi_integer,&
& glob_coo%ia,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info == 0) call mpi_allgatherv(loc_coo%ja,ndx,mpi_integer,&
& glob_coo%ja,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
if (keeploc_) then
call loca%mv_from(loc_coo)
else
call loc_coo%free()
end if
call glob_coo%set_nzeros(nzg)
if (present(dupl)) call glob_coo%set_dupl(dupl)
call globa%mv_from(glob_coo)
else
write(0,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dsp_allgather
base/modules/Makefile
+1
View File
@@ -53,6 +53,7 @@ psi_mod.o: psb_penv_mod.o psb_error_mod.o psb_desc_type.o psb_const_mod.o psi_se
psb_desc_type.o: psb_const_mod.o psb_error_mod.o psb_penv_mod.o psb_realloc_mod.o psb_hash_mod.o
psb_linmap_mod.o: psb_linmap_type_mod.o psb_mat_mod.o
psb_linmap_type_mod.o: psb_desc_type.o psb_error_mod.o psb_serial_mod.o psb_comm_mod.o psb_mat_mod.o
psb_comm_mod.o: psb_desc_type.o psb_mat_mod.o
psb_check_mod.o: psb_desc_type.o
psb_serial_mod.o: psb_mat_mod.o psb_string_mod.o psb_sort_mod.o psi_serial_mod.o
psb_sort_mod.o: psb_error_mod.o psb_realloc_mod.o psb_const_mod.o
base/modules/psb_comm_mod.f90
+11
View File
@@ -302,6 +302,17 @@ module psb_comm_mod
end interface
interface psb_gather
subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
use psb_descriptor_type
use psb_mat_mod
implicit none
type(psb_d_sparse_mat), intent(inout) :: loca
type(psb_d_sparse_mat), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer, intent(in), optional :: root,dupl
logical, intent(in), optional :: keepnum,keeploc
end subroutine psb_dsp_allgather
subroutine psb_igatherm(globx, locx, desc_a, info, root)
use psb_descriptor_type
integer, intent(in) :: locx(:,:)
base/modules/psb_d_mat_mod.f03
+143 -9
View File
@@ -52,18 +52,25 @@ module psb_d_mat_mod
procedure, pass(a) :: d_csgetrow
procedure, pass(a) :: d_csgetblk
generic, public :: csget => d_csgetptn, d_csgetrow, d_csgetblk
procedure, pass(a) :: csclip
procedure, pass(a) :: d_csclip
procedure, pass(a) :: d_b_csclip
generic, public :: csclip => d_b_csclip, d_csclip
procedure, pass(a) :: reall => reallocate_nz
procedure, pass(a) :: get_neigh
procedure, pass(a) :: d_cscnv
procedure, pass(a) :: d_cscnv_ip
generic, public :: cscnv => d_cscnv, d_cscnv_ip
procedure, pass(a) :: d_cscnv_base
generic, public :: cscnv => d_cscnv, d_cscnv_ip, d_cscnv_base
procedure, pass(a) :: reinit
procedure, pass(a) :: print => sparse_print
procedure, pass(a) :: d_mv_from
generic, public :: mv_from => d_mv_from
procedure, pass(a) :: d_mv_to
generic, public :: mv_to => d_mv_to
procedure, pass(a) :: d_cp_from
generic, public :: cp_from => d_cp_from
procedure, pass(a) :: d_cp_to
generic, public :: cp_to => d_cp_to
procedure, pass(a) :: d_transp_1mat
procedure, pass(a) :: d_transp_2mat
generic, public :: transp => d_transp_1mat, d_transp_2mat
@@ -92,7 +99,7 @@ module psb_d_mat_mod
& get_state, get_dupl, is_null, is_bld, is_upd, &
& is_asb, is_sorted, is_upper, is_lower, is_triangle, &
& is_unit, get_neigh, csall, csput, d_csgetrow,&
& d_csgetblk, csclip, d_cscnv, d_cscnv_ip, &
& d_csgetblk, d_csclip, d_b_csclip, d_cscnv, d_cscnv_ip, &
& reallocate_nz, free, trim, &
& sparse_print, reinit, &
& set_nrows, set_ncols, set_dupl, &
@@ -100,7 +107,7 @@ module psb_d_mat_mod
& set_upd, set_asb, set_sorted, &
& set_upper, set_lower, set_triangle, &
& set_unit, get_diag, get_nz_row, d_csgetptn, &
& d_mv_from, d_cp_from, &
& d_mv_from, d_mv_to, d_cp_from, d_cp_to,&
& d_transp_1mat, d_transp_2mat, &
& d_transc_1mat, d_transc_2mat
@@ -1262,7 +1269,7 @@ contains
subroutine csclip(a,b,info,&
subroutine d_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
@@ -1306,7 +1313,52 @@ contains
return
end if
end subroutine csclip
end subroutine d_csclip
subroutine d_b_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod
implicit none
class(psb_d_sparse_mat), intent(in) :: a
type(psb_d_coo_sparse_mat), intent(out) :: b
integer,intent(out) :: info
integer, intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
Integer :: err_act
character(len=20) :: name='csclip'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat), allocatable :: acoo
info = 0
call psb_erractionsave(err_act)
if (a%is_null()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == 0) call a%a%csclip(b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine d_b_csclip
@@ -1494,6 +1546,61 @@ contains
end subroutine d_cscnv_ip
subroutine d_cscnv_base(a,b,info,dupl)
use psb_error_mod
use psb_string_mod
implicit none
class(psb_d_sparse_mat), intent(in) :: a
class(psb_d_base_sparse_mat), intent(out) :: b
integer, intent(out) :: info
integer,optional, intent(in) :: dupl
type(psb_d_coo_sparse_mat) :: altmp
Integer :: err_act
character(len=20) :: name='cscnv'
logical, parameter :: debug=.false.
info = 0
call psb_erractionsave(err_act)
if (a%is_null()) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call a%a%cp_to_coo(altmp,info )
if ((info == 0).and.present(dupl)) then
call altmp%set_dupl(dupl)
end if
call altmp%fix(info)
if (info == 0) call altmp%trim()
if (info == 0) call altmp%set_asb()
if (info == 0) call b%mv_from_coo(altmp,info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name,a_err="mv_from")
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine d_cscnv_base
subroutine d_mv_from(a,b)
use psb_error_mod
use psb_string_mod
@@ -1538,6 +1645,33 @@ contains
end if
end subroutine d_cp_from
subroutine d_mv_to(a,b)
use psb_error_mod
use psb_string_mod
implicit none
class(psb_d_sparse_mat), intent(inout) :: a
class(psb_d_base_sparse_mat), intent(out) :: b
integer :: info
call b%mv_from_fmt(a%a,info)
return
end subroutine d_mv_to
subroutine d_cp_to(a,b)
use psb_error_mod
use psb_string_mod
implicit none
class(psb_d_sparse_mat), intent(in) :: a
class(psb_d_base_sparse_mat), intent(out) :: b
integer :: info
call b%cp_from_fmt(a%a,info)
return
end subroutine d_cp_to
subroutine d_sparse_mat_move(a,b,info)
use psb_error_mod
use psb_string_mod
@@ -1979,7 +2113,7 @@ contains
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='csnmi'
character(len=20) :: name='get_diag'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
@@ -2015,7 +2149,7 @@ contains
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='csnmi'
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
@@ -2051,7 +2185,7 @@ contains
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='csnmi'
character(len=20) :: name='scal'
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
base/modules/psb_serial_mod.f90
+54
View File
@@ -41,5 +41,59 @@ module psb_serial_mod
& psb_sct => psi_sct
use psb_mat_mod
interface psb_symbmm
subroutine psb_dsymbmm(a,b,c,info)
use psb_mat_mod
implicit none
type(psb_d_sparse_mat), intent(in) :: a,b
type(psb_d_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_dsymbmm
subroutine psb_dbase_symbmm(a,b,c,info)
use psb_mat_mod
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
end subroutine psb_dbase_symbmm
end interface
interface psb_numbmm
subroutine psb_dnumbmm(a,b,c)
use psb_mat_mod
implicit none
type(psb_d_sparse_mat), intent(in) :: a,b
type(psb_d_sparse_mat), intent(inout) :: c
end subroutine psb_dnumbmm
subroutine psb_dbase_numbmm(a,b,c)
use psb_mat_mod
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(inout) :: c
end subroutine psb_dbase_numbmm
end interface
interface psb_rwextd
subroutine psb_drwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
type(psb_d_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_d_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_drwextd
subroutine psb_dbase_rwextd(nr,a,info,b,rowscale)
use psb_mat_mod
implicit none
integer, intent(in) :: nr
class(psb_d_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
end subroutine psb_dbase_rwextd
end interface
end module psb_serial_mod
base/serial/Makefile
+2 -1
View File
@@ -1,7 +1,8 @@
include ../../Make.inc
FOBJS = psb_lsame.o
FOBJS = psb_lsame.o psb_dsymbmm.o psb_dnumbmm.o \
psb_drwextd.o
# psb_regen_mod.o psb_lsame.o psb_zspgetrow.o\
# psb_zcsmm.o psb_zcsmv.o psb_zspgtdiag.o psb_zspgtblk.o\
# psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\
base/serial/f77/Makefile
+1 -1
View File
@@ -4,7 +4,7 @@ include ../../../Make.inc
# The object files
#
FOBJS = daxpby.o saxpby.o \
caxpby.o zaxpby.o
caxpby.o zaxpby.o smmp.o
# clpupd.o ccsmm.o cswmm.o ccsnmi.o ccsrws.o\
base/serial/psb_dnumbmm.f90
+142 -41
View File
@@ -41,45 +41,146 @@
!
subroutine psb_dnumbmm(a,b,c)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dnumbmm
implicit none
implicit none
type(psb_dspmat_type) :: a,b,c
real(psb_dpk_), allocatable :: temp(:)
integer :: info
logical :: csra, csrb
allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
if (info /= 0) then
return
type(psb_d_sparse_mat), intent(in) :: a,b
type(psb_d_sparse_mat), intent(inout) :: c
integer :: info
integer :: err_act
character(len=*), parameter :: name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(aa=>c%a)
type is (psb_d_csr_sparse_mat)
call psb_numbmm(a%a,b%a,aa)
class default
info = 1121
call psb_errpush(info,name)
goto 9999
end select
call c%set_asb()
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dnumbmm
subroutine psb_dbase_numbmm(a,b,c)
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dbase_numbmm
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(inout) :: c
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
real(psb_dpk_), allocatable :: temp(:)
integer :: info
integer :: err_act
name='psb_numbmm'
call psb_erractionsave(err_act)
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(temp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
!
! Note: we still have to test about possible performance hits.
!
!
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
if (csra.and.csrb) then
call numbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
& b%ia2,b%ia1,0,b%aspk,&
& c%ia2,c%ia1,0,c%aspk,temp)
else
call inner_numbmm(a,b,c,temp,info)
if (info /= 0) then
write(0,*) 'Error ',info,' from inner numbmm'
end if
call psb_ensure_size(size(c%ja),c%val,info)
select type(a)
type is (psb_d_csr_sparse_mat)
select type(b)
type is (psb_d_csr_sparse_mat)
call csr_numbmm(a,b,c,temp,info)
class default
call gen_numbmm(a,b,c,temp,info)
end select
class default
call gen_numbmm(a,b,c,temp,info)
end select
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_sp_setifld(psb_spmat_asb_,psb_state_,c,info)
call c%set_asb()
deallocate(temp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
contains
subroutine inner_numbmm(a,b,c,temp,info)
type(psb_dspmat_type) :: a,b,c
subroutine csr_numbmm(a,b,c,temp,info)
type(psb_d_csr_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(inout) :: c
real(psb_dpk_) :: temp(:)
integer, intent(out) :: info
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
call numbmm(ma,na,nb,a%irp,a%ja,0,a%val,&
& b%irp,b%ja,0,b%val,&
& c%irp,c%ja,0,c%val,temp)
end subroutine csr_numbmm
subroutine gen_numbmm(a,b,c,temp,info)
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(inout) :: c
integer :: info
real(psb_dpk_) :: temp(:)
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
@@ -87,13 +188,14 @@ contains
integer :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
real(psb_dpk_) :: ajj
n = a%m
m = a%k
l = b%k
n = a%get_nrows()
m = a%get_ncols()
l = b%get_ncols()
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
if (info /= 0) then
if (info /= 0) then
info = 4000
return
endif
@@ -105,8 +207,7 @@ contains
minmn = min(m,n)
do i = 1,n
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
call a%csget(i,i,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
ajj = aval(jj)
@@ -116,7 +217,7 @@ contains
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
@@ -127,19 +228,19 @@ contains
endif
enddo
end do
do j = c%ia2(i),c%ia2(i+1)-1
if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
do j = c%irp(i),c%irp(i+1)-1
if((c%ja(j)<1).or. (c%ja(j) > maxlmn)) then
write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
info = 3
return
else
c%aspk(j) = temp(c%ia1(j))
temp(c%ia1(j)) = dzero
c%val(j) = temp(c%ja(j))
temp(c%ja(j)) = dzero
endif
end do
end do
end subroutine gen_numbmm
end subroutine inner_numbmm
end subroutine psb_dnumbmm
end subroutine psb_dbase_numbmm
base/serial/psb_drwextd.f90
+182 -100
View File
@@ -39,127 +39,54 @@
!
!
subroutine psb_drwextd(nr,a,info,b,rowscale)
use psb_spmat_type
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_drwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
type(psb_dspmat_type), intent(inout) :: a
integer,intent(out) :: info
type(psb_dspmat_type), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer, intent(in) :: nr
type(psb_d_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
type(psb_d_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb
character(len=20) :: name, ch_err
type(psb_d_coo_sparse_mat) :: actmp
logical rowscale_
name='psb_drwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
if (nr > a%m) then
if (psb_toupper(a%fida) == 'CSR') then
call psb_ensure_size(nr+1,a%ia2,info)
if (nr > a%get_nrows()) then
select type(aa=> a%a)
type is (psb_d_csr_sparse_mat)
if (present(b)) then
nzb = psb_sp_get_nnzeros(b)
call psb_ensure_size(size(a%ia1)+nzb,a%ia1,info)
call psb_ensure_size(size(a%aspk)+nzb,a%aspk,info)
if (psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
a%ia2(a%m+i+1) = a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i)
ja = a%ia2(a%m+i)
jb = b%ia2(i)
do
if (jb >= b%ia2(i+1)) exit
a%aspk(ja) = b%aspk(jb)
a%ia1(ja) = b%ia1(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-a%m
a%ia2(a%m+i+1) = a%ia2(a%m+i)
end do
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
call psb_rwextd(nr,aa,info,b%a,rowscale)
else
do i=a%m+2,nr+1
a%ia2(i) = a%ia2(i-1)
end do
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
a%m = nr
a%k = max(a%k,b%k)
else if (psb_toupper(a%fida) == 'COO') then
type is (psb_d_coo_sparse_mat)
if (present(b)) then
nza = psb_sp_get_nnzeros(a)
nzb = psb_sp_get_nnzeros(b)
call psb_sp_reall(a,nza+nzb,info)
if (psb_toupper(b%fida)=='COO') then
if (rowscale_) then
do j=1,nzb
if ((a%m + b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = a%m + b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
end if
enddo
else
do j=1,nzb
if ((b%ia1(j)) <= nr) then
nza = nza + 1
a%ia1(nza) = b%ia1(j)
a%ia2(nza) = b%ia2(j)
a%aspk(nza) = b%aspk(j)
endif
enddo
endif
a%infoa(psb_nnz_) = nza
else if(psb_toupper(b%fida)=='CSR') then
do i=1, min(nr-a%m,b%m)
do
jb = b%ia2(i)
if (jb >= b%ia2(i+1)) exit
nza = nza + 1
a%aspk(nza) = b%aspk(jb)
a%ia1(nza) = a%m + i
a%ia2(nza) = b%ia1(jb)
jb = jb + 1
end do
end do
a%infoa(psb_nnz_) = nza
call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
else
call psb_rwextd(nr,aa,info,rowscale=rowscale)
end if
class default
call aa%mv_to_coo(actmp,info)
if (info == 0) then
if (present(b)) then
call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
else
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
endif
endif
a%m = nr
a%k = max(a%k,b%k)
else if (a%fida == 'JAD') then
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
else
info=136
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_rwextd(nr,actmp,info,rowscale=rowscale)
end if
end if
if (info == 0) call aa%mv_from_coo(actmp,info)
end select
end if
if (info /= 0) goto 9999
call psb_erractionrestore(err_act)
return
@@ -173,3 +100,158 @@ subroutine psb_drwextd(nr,a,info,b,rowscale)
return
end subroutine psb_drwextd
subroutine psb_dbase_rwextd(nr,a,info,b,rowscale)
use psb_error_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dbase_rwextd
implicit none
! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
integer, intent(in) :: nr
class(psb_d_base_sparse_mat), intent(inout) :: a
integer,intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: b
logical,intent(in), optional :: rowscale
integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
character(len=20) :: name, ch_err
logical rowscale_
name='psb_dbase_rwextd'
info = 0
call psb_erractionsave(err_act)
if (present(rowscale)) then
rowscale_ = rowscale
else
rowscale_ = .true.
end if
ma = a%get_nrows()
na = a%get_ncols()
select type(a)
type is (psb_d_csr_sparse_mat)
call psb_ensure_size(nr+1,a%irp,info)
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
select type (b)
type is (psb_d_csr_sparse_mat)
call psb_ensure_size(size(a%ja)+nzb,a%ja,info)
call psb_ensure_size(size(a%val)+nzb,a%val,info)
do i=1, min(nr-ma,mb)
a%irp(ma+i+1) = a%irp(ma+i) + b%irp(i+1) - b%irp(i)
ja = a%irp(ma+i)
jb = b%irp(i)
do
if (jb >= b%irp(i+1)) exit
a%val(ja) = b%val(jb)
a%ja(ja) = b%ja(jb)
ja = ja + 1
jb = jb + 1
end do
end do
do j=i,nr-ma
a%irp(ma+i+1) = a%irp(ma+i)
end do
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
else
do i=ma+2,nr+1
a%irp(i) = a%irp(i-1)
end do
end if
call a%set_nrows(nr)
type is (psb_d_coo_sparse_mat)
nza = a%get_nzeros()
if (present(b)) then
mb = b%get_nrows()
nb = b%get_ncols()
nzb = b%get_nzeros()
call a%reallocate(nza+nzb)
select type(b)
type is (psb_d_coo_sparse_mat)
if (rowscale_) then
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = ma + b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
else
do j=1,nzb
if ((ma + b%ia(j)) <= nr) then
nza = nza + 1
a%ia(nza) = b%ia(j)
a%ja(nza) = b%ja(j)
a%val(nza) = b%val(j)
end if
enddo
endif
call a%set_nzeros(nza)
type is (psb_d_csr_sparse_mat)
do i=1, min(nr-ma,mb)
do
jb = b%irp(i)
if (jb >= b%irp(i+1)) exit
nza = nza + 1
a%val(nza) = b%val(jb)
a%ia(nza) = ma + i
a%ja(nza) = b%ja(jb)
jb = jb + 1
end do
end do
call a%set_nzeros(nza)
class default
write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
end select
call a%set_ncols(max(na,nb))
endif
call a%set_nrows(nr)
class default
info = 135
ch_err=a%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dbase_rwextd
base/serial/psb_dsymbmm.f90
+173 -95
View File
@@ -40,62 +40,107 @@
!
subroutine psb_dsymbmm(a,b,c,info)
use psb_spmat_type
use psb_mat_mod
use psb_string_mod
use psb_serial_mod, psb_protect_name => psb_dsymbmm
implicit none
type(psb_dspmat_type) :: a,b,c
integer, allocatable :: itemp(:)
integer :: nze,info
type(psb_d_sparse_mat), intent(in) :: a,b
type(psb_d_sparse_mat), intent(out) :: c
integer, intent(out) :: info
type(psb_d_csr_sparse_mat), allocatable :: ccsr
integer :: err_act
character(len=*), parameter :: name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
if ((a%is_null()) .or.(b%is_null())) then
info = 1121
call psb_errpush(info,name)
goto 9999
endif
allocate(ccsr, stat=info)
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call psb_symbmm(a%a,b%a,ccsr,info)
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call move_alloc(ccsr,c%a)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dsymbmm
subroutine psb_dbase_symbmm(a,b,c,info)
use psb_mat_mod
use psb_serial_mod, psb_protect_name => psb_dbase_symbmm
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(out) :: c
integer, intent(out) :: info
integer, allocatable :: itemp(:)
integer :: nze, ma,na,mb,nb
character(len=20) :: name
integer :: err_act
logical :: csra, csrb
name='psb_symbmm'
call psb_erractionsave(err_act)
info = 0
csra = (psb_toupper(a%fida(1:3))=='CSR')
csrb = (psb_toupper(b%fida(1:3))=='CSR')
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
if ( mb /= na ) then
write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
endif
allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info)
allocate(itemp(max(ma,na,mb,nb)),stat=info)
if (info /= 0) then
return
info = 4000
call psb_Errpush(info,name)
goto 9999
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
!
if (csra.and.csrb) then
call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
& b%ia2,b%ia1,0,&
& c%ia2,c%ia1,0,itemp)
else
call inner_symbmm(a,b,c,itemp,info)
endif
call psb_realloc(size(c%ia1),c%aspk,info)
select type(a)
type is (psb_d_csr_sparse_mat)
select type(b)
type is (psb_d_csr_sparse_mat)
call csr_symbmm(a,b,c,itemp,info)
class default
call gen_symbmm(a,b,c,itemp,info)
end select
class default
call gen_symbmm(a,b,c,itemp,info)
end select
c%pl(1) = 0
c%pr(1) = 0
c%m=a%m
c%k=b%k
c%fida='CSR'
c%descra='GUN'
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(size(c%ja),c%val,info)
deallocate(itemp)
call psb_erractionrestore(err_act)
return
@@ -109,86 +154,119 @@ subroutine psb_dsymbmm(a,b,c,info)
contains
subroutine inner_symbmm(a,b,c,index,info)
type(psb_dspmat_type) :: a,b,c
subroutine csr_symbmm(a,b,c,itemp,info)
type(psb_d_csr_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(out) :: c
integer :: itemp(:)
integer, intent(out) :: info
interface
subroutine symbmm (n, m, l, ia, ja, diaga, &
& ib, jb, diagb, ic, jc, diagc, index)
integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& diagc, index(*)
integer, allocatable :: ic(:),jc(:)
end subroutine symbmm
end interface
integer :: nze, ma,na,mb,nb
info = 0
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
call symbmm(ma,na,nb,a%irp,a%ja,0,&
& b%irp,b%ja,0,&
& c%irp,c%ja,0,itemp)
end subroutine csr_symbmm
subroutine gen_symbmm(a,b,c,index,info)
class(psb_d_base_sparse_mat), intent(in) :: a,b
type(psb_d_csr_sparse_mat), intent(out) :: c
integer :: index(:),info
integer, allocatable :: iarw(:), iacl(:),ibrw(:),ibcl(:)
real(psb_dpk_), allocatable :: aval(:),bval(:)
integer :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
integer :: nze, ma,na,mb,nb
ma = a%get_nrows()
na = a%get_ncols()
mb = b%get_nrows()
nb = b%get_ncols()
n = a%m
m = a%k
l = b%k
nze = max(ma+1,2*ma)
call c%allocate(ma,nb,nze)
n = ma
m = na
l = nb
maxlmn = max(l,m,n)
allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
& aval(maxlmn),bval(maxlmn), stat=info)
if (info /= 0) then
info = 4000
return
endif
if (size(c%ia2) < n+1) then
call psb_realloc(n+1,c%ia2,info)
endif
do i=1,maxlmn
index(i)=0
end do
c%ia2(1)=1
minlm = min(l,m)
minmn = min(m,n)
c%irp(1)=1
minlm = min(l,m)
minmn = min(m,n)
main: do i=1,n
istart=-1
length=0
call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
main: do i=1,n
istart=-1
length=0
call a%csget(i,i,nazr,iarw,iacl,aval,info)
do jj=1, nazr
j=iacl(jj)
if ((j<1).or.(j>m)) then
write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
info = 1
return
endif
call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
endif
call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
do k=1,nbzr
if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then
write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
info=2
return
else
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
end do
c%ia2(i+1)=c%ia2(i)+length
if (c%ia2(i+1) > size(c%ia1)) then
if (n > (2*i)) then
nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i))
else
nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ia1,info)
end if
do j= c%ia2(i),c%ia2(i+1)-1
c%ia1(j)=istart
istart=index(istart)
index(c%ia1(j))=0
if(index(ibcl(k)) == 0) then
index(ibcl(k))=istart
istart=ibcl(k)
length=length+1
endif
endif
end do
call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
index(i) = 0
end do main
end do
end subroutine inner_symbmm
c%irp(i+1)=c%irp(i)+length
end subroutine psb_dsymbmm
if (c%irp(i+1) > size(c%ja)) then
if (n > (2*i)) then
nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
else
nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,c%ja,info)
end if
do j= c%irp(i),c%irp(i+1)-1
c%ja(j)=istart
istart=index(istart)
index(c%ja(j))=0
end do
call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
index(i) = 0
end do main
end subroutine gen_symbmm
end subroutine psb_dbase_symbmm
configure
Vendored
+11 -13
View File
@@ -771,7 +771,7 @@ SHELL'
ac_subst_files=''
ac_user_opts='
enable_option_checking
with_serial_mpi
enable_serial
with_ccopt
with_fcopt
with_f90copt
@@ -1429,14 +1429,14 @@ Optional Features:
--disable-option-checking ignore unrecognized --enable/--with options
--disable-FEATURE do not include FEATURE (same as --enable-FEATURE=no)
--enable-FEATURE[=ARG] include FEATURE [ARG=yes]
--enable-serial Specify whether to enable a fake mpi library to run
in serial mode.
--disable-dependency-tracking speeds up one-time build
--enable-dependency-tracking do not reject slow dependency extractors
Optional Packages:
--with-PACKAGE[=ARG] use PACKAGE [ARG=yes]
--without-PACKAGE do not use PACKAGE (same as --with-PACKAGE=no)
--with-serial-mpi Specify whether to enable a fake mpi library to run
in serial mode. --with-serial-mpi={yes|no}
--with-ccopt additional CCOPT flags to be added: will prepend
to CCOPT
--with-fcopt additional FCOPT flags to be added: will prepend
@@ -3394,14 +3394,11 @@ fi
{ $as_echo "$as_me:$LINENO: checking whether we want serial (fake) mpi" >&5
$as_echo_n "checking whether we want serial (fake) mpi... " >&6; }
# Check whether --enable-serial was given.
if test "${enable_serial+set}" = set; then
enableval=$enable_serial;
pac_cv_serial_mpi="yes";
# Check whether --with-serial-mpi was given.
if test "${with_serial_mpi+set}" = set; then
withval=$with_serial_mpi;
pac_cv_serial_mpi="${withval}";
else
pac_cv_serial_mpi="no";
fi
@@ -3409,7 +3406,8 @@ if test x"$pac_cv_serial_mpi" == x"yes" ; then
{ $as_echo "$as_me:$LINENO: result: yes." >&5
$as_echo "yes." >&6; }
else
{ $as_echo "$as_me:$LINENO: result: no." >&5
pac_cv_serial_mpi="no";
{ $as_echo "$as_me:$LINENO: result: no." >&5
$as_echo "no." >&6; }
fi
@@ -7119,7 +7117,7 @@ if test "X$FCOPT" == "X" ; then
fi
if test "X$psblas_cv_fc" == X"nag" ; then
# Add needed options
FCOPT="$FCOPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
FCOPT="$FCOPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
fi
FFLAGS="${FCOPT}"
@@ -7152,7 +7150,7 @@ else
fi
if test "X$psblas_cv_fc" == X"nag" ; then
# Add needed options
F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
EXTRA_OPT="-mismatch_all"
F03COPT="${F90COPT}"
else
configure.ac
+2 -2
View File
@@ -399,7 +399,7 @@ if test "X$FCOPT" == "X" ; then
fi
if test "X$psblas_cv_fc" == X"nag" ; then
# Add needed options
FCOPT="$FCOPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
FCOPT="$FCOPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
fi
FFLAGS="${FCOPT}"
@@ -432,7 +432,7 @@ else
fi
if test "X$psblas_cv_fc" == X"nag" ; then
# Add needed options
F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
EXTRA_OPT="-mismatch_all"
F03COPT="${F90COPT}"
else