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

*** empty log message ***

Browse Source
This commit is contained in:
Alfredo Buttari
2005-10-10 08:33:44 +00:00
parent 3ca7ba61e5
commit d2afa63600
21 changed files with 497 additions and 1192 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/comm/psb_dgather.f90
+4 -4
View File
@@ -243,7 +243,7 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,&
root = -1
end if
if (root==-1) then
iiroot=0
root=0
endif
jglobx=1
@@ -268,10 +268,10 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,&
k = 1
if (myrow == iiroot) then
if (myrow == root) then
call igebs2d(icontxt, 'all', ' ', 1, 1, k, 1)
else
call igebr2d(icontxt, 'all', ' ', 1, 1, k, 1, iiroot, 0)
call igebr2d(icontxt, 'all', ' ', 1, 1, k, 1, root, 0)
end if
! there should be a global check on k here!!!
@@ -298,7 +298,7 @@ subroutine psb_dgatherv(globx, locx, desc_a, info, iroot,&
globx(idx) = locx(i)
end do
! adjust overlapped elements
i=0
i=1
do while (desc_a%ovrlap_elem(i).ne.-1)
idx=desc_a%ovrlap_elem(i+psb_ovrlp_elem_)
idx=desc_a%loc_to_glob(idx)
src/comm/psb_dhalo.f90
+6
View File
@@ -162,6 +162,9 @@ subroutine psb_dhalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
@@ -314,6 +317,9 @@ subroutine psb_dhalov(x,desc_a,info,alpha,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
src/comm/psb_dovrl.f90
+5 -3
View File
@@ -173,7 +173,8 @@ subroutine psb_dovrlm(x,desc_a,info,jx,ik,work,choice,update_type)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
@@ -349,8 +350,9 @@ subroutine psb_dovrlv(x,desc_a,info,work,choice,update_type)
goto 9999
end select
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
src/comm/psb_ihalo.f90
+6
View File
@@ -158,6 +158,9 @@ subroutine psb_ihalom(x,desc_a,info,alpha,jx,ik,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
@@ -311,6 +314,9 @@ subroutine psb_ihalov(x,desc_a,info,alpha,work,tran,mode)
goto 9999
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
src/modules/Makefile
+1 -1
View File
@@ -6,7 +6,7 @@ MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \
psb_prec_type.o psb_error_mod.o psb_prec_mod.o \
psb_methd_mod.o psb_const_mod.o \
psb_comm_mod.o psb_psblas_mod.o psi_mod.o \
psb_sparse_mod.o psb_check_mod.o
psb_sparse_mod.o psb_check_mod.o psb_all_mod.o
OBJS = error.o
src/modules/psb_methd_mod.f90
+1 -1
View File
@@ -19,7 +19,7 @@ Module psb_methd_mod
end subroutine psb_dcg
end interface
interface spb_bicg
interface psb_bicg
subroutine psb_dbicg(a,prec,b,x,eps,&
& desc_a,info,itmax,iter,err,itrace,istop)
use psb_serial_mod
src/modules/psb_psblas_mod.f90
+2 -7
View File
@@ -74,9 +74,6 @@ module psb_psblas_mod
type(psb_desc_type), intent (in) :: desc_a
integer, intent(out) :: info
end function psb_damaxv
end interface
interface psb_amaxs
subroutine psb_damaxvs(res,x,desc_a,info)
use psb_descriptor_type
real(kind(1.d0)), intent (out) :: res
@@ -84,12 +81,13 @@ module psb_psblas_mod
type(psb_desc_type), intent (in) :: desc_a
integer, intent(out) :: info
end subroutine psb_damaxvs
subroutine psb_dmamax(res,x,desc_a,info)
subroutine psb_dmamax(res,x,desc_a,info,jx)
use psb_descriptor_type
real(kind(1.d0)), intent (out) :: res(:)
real(kind(1.d0)), intent (in) :: x(:,:)
type(psb_desc_type), intent (in) :: desc_a
integer, intent(out) :: info
integer, optional :: jx
end subroutine psb_dmamax
end interface
@@ -145,9 +143,6 @@ module psb_psblas_mod
type(psb_desc_type), intent (in) :: desc_a
integer, intent(out) :: info
end function psb_dnrm2v
end interface
interface psb_nrm2s
subroutine psb_dnrm2vs(res,x,desc_a,info)
use psb_descriptor_type
real(kind(1.d0)), intent (out) :: res
src/modules/psb_tools_mod.f90
+5 -5
View File
@@ -391,11 +391,11 @@ Module psb_tools_mod
type(psb_dspmat_type), intent(inout) ::a
integer, intent(out) :: info
end subroutine psb_dspfree
subroutine psb_dspfrees(a,info)
use psb_spmat_type
type(psb_dspmat_type), intent(inout) ::a
integer, intent(out) :: info
end subroutine psb_dspfrees
!!$ subroutine psb_dspfrees(a,info)
!!$ use psb_spmat_type
!!$ type(psb_dspmat_type), intent(inout) ::a
!!$ integer, intent(out) :: info
!!$ end subroutine psb_dspfrees
end interface
src/prec/psb_dprecfree.f90
+1 -1
View File
@@ -50,7 +50,7 @@ subroutine psb_dprecfree(p,info)
if (associated(p%baseprecv)) then
do i=1,size(p%baseprecv)
call psb_baseprecfree(p%baseprecv(i),info)
call psb_base_precfree(p%baseprecv(i),info)
end do
deallocate(p%baseprecv)
p%baseprecv => null()
src/psblas/psb_damax.f90
+13 -14
View File
@@ -217,7 +217,7 @@ end function psb_damaxv
! info - integer. Eventually returns an error code.
! jx - integer(optional). The column offset.
!
subroutine psb_damaxvs (res,x,desc_a, info, jx)
subroutine psb_damaxvs (res,x,desc_a, info)
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
@@ -227,14 +227,12 @@ subroutine psb_damaxvs (res,x,desc_a, info, jx)
real(kind(1.d0)), intent(in) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer, optional, intent(in) :: jx
real(kind(1.D0)), intent(out) :: res
! locals
integer :: int_err(5), icontxt, nprow, npcol, myrow, mycol,&
& err_act, n, iix, jjx, temp(2), ix, ijx, m, imax, idamax
real(kind(1.d0)) :: locmax(2), amax
real(kind(1.d0)),pointer :: tmpx(:)
character(len=20) :: name, ch_err
name='psb_damaxvs'
@@ -258,13 +256,9 @@ subroutine psb_damaxvs (res,x,desc_a, info, jx)
call psb_errpush(info,name)
goto 9999
endif
ix = 1
if (present(jx)) then
ijx = jx
else
ijx = 1
endif
ijx=1
m = desc_a%matrix_data(psb_m_)
@@ -285,7 +279,7 @@ subroutine psb_damaxvs (res,x,desc_a, info, jx)
! compute local max
if ((desc_a%matrix_data(psb_n_row_).gt.0).and.(m.ne.0)) then
imax=idamax(desc_a%matrix_data(psb_n_row_)-iix+1,x(iix),1)
amax=abs(tmpx(iix+imax-1))
amax=abs(x(iix+imax-1))
end if
! compute global max
@@ -321,7 +315,7 @@ end subroutine psb_damaxvs
! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code.
!
subroutine psb_dmamaxs (res,x,desc_a, info)
subroutine psb_dmamaxs (res,x,desc_a, info,jx)
use psb_serial_mod
use psb_descriptor_type
use psb_check_mod
@@ -331,11 +325,12 @@ subroutine psb_dmamaxs (res,x,desc_a, info)
real(kind(1.d0)), intent(in) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info
integer, optional, intent(in) :: jx
real(kind(1.d0)), intent(out) :: res(:)
! locals
integer :: int_err(5), icontxt, nprow, npcol, myrow, mycol,&
& err_act, n, iix, jjx, ix, jx, temp(2), ijx, m, imax, i, k, idamax
& err_act, n, iix, jjx, ix, temp(2), ijx, m, imax, i, k, idamax
real(kind(1.d0)) :: locmax(2), amax
real(kind(1.d0)),pointer :: tmpx(:)
character(len=20) :: name, ch_err
@@ -363,12 +358,16 @@ subroutine psb_dmamaxs (res,x,desc_a, info)
endif
ix = 1
jx = 1
if (present(jx)) then
ijx = jx
else
ijx = 1
endif
m = desc_a%matrix_data(psb_m_)
k = min(size(x,2),size(res,1))
call psb_chkvect(m,1,size(x,1),ix,jx,desc_a%matrix_data,info,iix,jjx)
call psb_chkvect(m,1,size(x,1),ix,ijx,desc_a%matrix_data,info,iix,jjx)
if(info.ne.0) then
info=4010
ch_err='psb_chkvect'
src/psblas/psb_dspmm.f90
+6 -1
View File
@@ -322,6 +322,7 @@ subroutine psb_dspmm(alpha,a,x,beta,y,desc_a,info,&
end if
if(.not.present(work)) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
@@ -455,6 +456,7 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
liwork= 2*ncol
if (a%pr(1) /= 0) liwork = liwork + n * ik
if (a%pl(1) /= 0) liwork = liwork + m * ik
! write(0,*)'---->>>',work(1)
if (present(work)) then
if(size(work).ge.liwork) then
iwork => work
@@ -593,7 +595,10 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
end if
if(.not.present(work)) deallocate(iwork)
if(.not.present(work)) then
deallocate(iwork)
end if
nullify(iwork)
call psb_erractionrestore(err_act)
return
src/tools/psb_dspfree.f90
+75 -75
View File
@@ -95,78 +95,78 @@ end subroutine psb_dspfree
subroutine psb_dspfrees(a, info)
!...free sparse matrix structure...
use psb_descriptor_type
use psb_spmat_type
use psb_serial_mod
use psb_const_mod
use psb_error_mod
implicit none
!....parameters...
type(psb_dspmat_type), intent(inout) ::a
integer, intent(out) :: info
!...locals....
integer :: int_err(5)
integer :: temp(1)
real(kind(1.d0)) :: real_err(5)
integer :: icontxt,nprow,npcol,me,mypcol,err, err_act
integer,parameter :: ione=1
character(len=20) :: name, ch_err
info=0
name = 'psb_dspfrees'
call psb_erractionsave(err_act)
!...deallocate a....
if ((info.eq.0).and.(.not.associated(a%pr))) info=2951
if (info.eq.0) then
!deallocate pr field
deallocate(a%pr,stat=info)
if (info.ne.0) info=2045
end if
if ((info.eq.0).and.(.not.associated(a%pl))) info=2952
!deallocate pl field
if (info.eq.0) then
deallocate(a%pl,stat=info)
if (info.ne.0) info=2046
end if
if ((info.eq.0).and.(.not.associated(a%ia2))) info=2953
if (info.eq.0) then
!deallocate ia2 field
deallocate(a%ia2,stat=info)
if (info.ne.0) info=2047
end if
if ((info.eq.0).and.(.not.associated(a%ia1))) info=2954
if (info.eq.0) then
!deallocate ia1 field
deallocate(a%ia1,stat=info)
if (info.ne.0) info=2048
endif
if ((info.eq.0).and.(.not.associated(a%aspk))) info=2955
if (info.eq.0) then
!deallocate aspk field
deallocate(a%aspk,stat=info)
if (info.ne.0) info=2049
endif
if (info.eq.0) call psb_nullify_sp(a)
if(info.ne.0) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error()
return
end if
return
end subroutine psb_dspfrees
!!$subroutine psb_dspfrees(a, info)
!!$ !...free sparse matrix structure...
!!$ use psb_descriptor_type
!!$ use psb_spmat_type
!!$ use psb_serial_mod
!!$ use psb_const_mod
!!$ use psb_error_mod
!!$ implicit none
!!$
!!$ !....parameters...
!!$ type(psb_dspmat_type), intent(inout) ::a
!!$ integer, intent(out) :: info
!!$ !...locals....
!!$ integer :: int_err(5)
!!$ integer :: temp(1)
!!$ real(kind(1.d0)) :: real_err(5)
!!$ integer :: icontxt,nprow,npcol,me,mypcol,err, err_act
!!$ integer,parameter :: ione=1
!!$ character(len=20) :: name, ch_err
!!$
!!$ info=0
!!$ name = 'psb_dspfrees'
!!$ call psb_erractionsave(err_act)
!!$
!!$ !...deallocate a....
!!$
!!$ if ((info.eq.0).and.(.not.associated(a%pr))) info=2951
!!$ if (info.eq.0) then
!!$ !deallocate pr field
!!$ deallocate(a%pr,stat=info)
!!$ if (info.ne.0) info=2045
!!$ end if
!!$ if ((info.eq.0).and.(.not.associated(a%pl))) info=2952
!!$ !deallocate pl field
!!$ if (info.eq.0) then
!!$ deallocate(a%pl,stat=info)
!!$ if (info.ne.0) info=2046
!!$ end if
!!$ if ((info.eq.0).and.(.not.associated(a%ia2))) info=2953
!!$ if (info.eq.0) then
!!$ !deallocate ia2 field
!!$ deallocate(a%ia2,stat=info)
!!$ if (info.ne.0) info=2047
!!$ end if
!!$ if ((info.eq.0).and.(.not.associated(a%ia1))) info=2954
!!$ if (info.eq.0) then
!!$ !deallocate ia1 field
!!$ deallocate(a%ia1,stat=info)
!!$ if (info.ne.0) info=2048
!!$ endif
!!$ if ((info.eq.0).and.(.not.associated(a%aspk))) info=2955
!!$ if (info.eq.0) then
!!$ !deallocate aspk field
!!$ deallocate(a%aspk,stat=info)
!!$ if (info.ne.0) info=2049
!!$ endif
!!$ if (info.eq.0) call psb_nullify_sp(a)
!!$
!!$ if(info.ne.0) then
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_erractionrestore(err_act)
!!$ return
!!$
!!$9999 continue
!!$ call psb_erractionrestore(err_act)
!!$ if (err_act.eq.act_abort) then
!!$ call psb_error()
!!$ return
!!$ end if
!!$ return
!!$
!!$end subroutine psb_dspfrees
test/Fileread/Makefile
+4 -82
View File
@@ -2,15 +2,8 @@ include ../../Make.inc
#
# Libraries used
#
LIBDIR=../../LIB/
LIBDIR=../../lib/
PSBLAS_LIB= -L$(LIBDIR) -lpsblas
SPARKER_LIB= -L$(LIBDIR) -lsparker
BLAS90LIB=-L$(LIBDIR) -lpsblas90 $(SLU)
#METHD90LIB=-L$(LIBDIR) -lmethd90
#TOOLS90LIB=-L$(LIBDIR) -ltools90
#PREC90LIB=-L$(LIBDIR) -lprec90
#
# We are using the public domain tool METIS from U. Minnesota. To get it
@@ -18,93 +11,22 @@ BLAS90LIB=-L$(LIBDIR) -lpsblas90 $(SLU)
#
METIS_LIB = -L$(HOME)/NUMERICAL/metis-4.0 -lmetis
#
CCOPT= -g
INCDIRS=-I$(LIBDIR)
TMMOBJS=partgraph.o part_block.o read_mat.o getp.o \
mmio.o mat_dist.o testmm.o
DFOBJS=partgraph.o part_block.o read_mat.o getp.o \
mmio.o mat_dist.o df_sample.o lowerc.o part_blk2.o
DAOBJS=partgraph.o part_block.o getp.o \
mmio.o mat_dist.o read_mat.o d_aggr.o part_blk2.o lowerc.o
DFLOBJS=partgraph.o part_block.o read_mat.o getp.o \
mmio.o mat_dist.o df_samplelog.o
ZFOBJS=partgraph.o part_block.o read_mat.o getp.o \
mmio.o mat_dist.o zf_sample.o
DFMOBJS=partgraph.o part_block.o mmio.o read_mat.o \
mat_dist.o df_samplem.o part_blk2.o lowerc.o
DFBOBJS=partgraph.o part_block.o mmio.o read_mat.o \
mat_dist.o comm_info.o df_bench.o part_blk2.o lowerc.o
DFCOBJS=partgraph.o part_block.o mmio.o read_mat.o \
mat_dist.o comm_info.o df_comm.o part_blk2.o lowerc.o
EXEDIR=./RUNS
all: df_sample zf_sample df_samplelog testmm df_samplem
all: df_sample
read_mat.o: mmio.o
df_sample: $(DFOBJS)
$(F90LINK) $(LINKOPT) $(DFOBJS) -o df_sample\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
$(PSBLAS_LIB) $(METIS_LIB) $(BLACS)
/bin/mv df_sample $(EXEDIR)
d_aggr: $(DAOBJS)
$(F90LINK) $(LINKOPT) $(DAOBJS) -o d_aggr\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
/bin/mv d_aggr $(EXEDIR)
df_samplem: $(DFMOBJS)
$(F90LINK) $(LINKOPT) $(DFMOBJS) -o df_samplem\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
/bin/mv df_samplem $(EXEDIR)
df_bench: $(DFBOBJS)
$(F90LINK) $(LINKOPT) $(DFBOBJS) -o df_bench\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
/bin/mv df_bench $(EXEDIR)
df_comm: $(DFCOBJS)
$(F90LINK) $(LINKOPT) $(DFCOBJS) -o df_comm\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
/bin/mv df_comm $(EXEDIR)
testmm: $(TMMOBJS)
$(F90LINK) $(LINKOPT) $(TMMOBJS) -o testmm\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(TOOLS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS)
/bin/mv testmm $(EXEDIR)
df_samplelog: $(DFLOBJS)
$(F90LINK) $(LINKOPT) $(DFLOBJS) -o df_samplelog\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(SPARKER_LIB) $(PREC90LIB) $(BLAS90LIB) $(PSBLAS_LIB) \
$(BLAS) $(SPARKER_LIB) $(BLACS) $(BLAS) -llmpe -lmpe
/bin/mv df_samplelog $(EXEDIR)
zf_sample: $(ZFOBJS)
$(F90LINK) $(LINKOPT) $(ZFOBJS) -o zf_sample\
$(METHD90LIB) $(TOOLS90LIB) $(BLAS90LIB) $(METIS_LIB)\
$(PSBLAS_LIB) $(SPARKER_LIB) $(BLAS)\
$(BLACS)
/bin/mv zf_sample $(EXEDIR)
aggr.o: mmio.o
aggr: aggr.o mmio.o
@@ -127,7 +49,7 @@ clean:
*$(.mod) $(EXEDIR)/df_sample $(EXEDIR)/zf_sample $(EXEDIR)/df_comm $(EXEDIR)/df_samplelog $(EXEDIR)/df_bench
lib:
(cd ../../; make lib)
(cd ../../; make library)
verycleanlib:
(cd ../../; make veryclean)
test/Fileread/RUNS/Makefile
+3 -8
View File
@@ -4,12 +4,7 @@ clean:
(cd ..; $(MAKE) clean)
verycleanlib:
(cd ..; $(MAKE) verycleanlib)
testmm:
(cd ..; $(MAKE) testmm)
df_sample:
(cd ..; $(MAKE) clean clean df_sample)
df_samplem:
(cd ..; $(MAKE) clean clean df_samplem)
zf_sample:
(cd ..; $(MAKE) clean zf_sample)
.PHONY: df_sample zf_sample testmm df_samplem
(cd ..; $(MAKE) df_sample)
.PHONY: df_sample
test/Fileread/RUNS/rtst.inp
+1 -1
View File
@@ -1,5 +1,5 @@
11 Number of inputs
kivap007.mtx This (and others) from: http://math.nist.gov/MatrixMarket/
kivap001.mtx This (and others) from: http://math.nist.gov/MatrixMarket/
NONE
BICGSTAB
ILU !!!! Actually, it's IPREC below. Should take this line out.
test/Fileread/df_sample.f90
+252 -246
View File
@@ -1,211 +1,211 @@
PROGRAM DF_SAMPLE
USE F90SPARSE
USE MAT_DIST
USE READ_MAT
USE PARTGRAPH
USE GETP
IMPLICIT NONE
program df_sample
use psb_all_mod
use mat_dist
use read_mat
use partgraph
use getp
implicit none
! Input parameters
CHARACTER*20 :: CMETHD, PREC, MTRX_FILE, RHS_FILE
CHARACTER*80 :: CHARBUF
! input parameters
character*20 :: cmethd, prec, mtrx_file, rhs_file
character*80 :: charbuf
DOUBLE PRECISION DDOT
EXTERNAL DDOT
INTERFACE
interface
! .....user passed subroutine.....
SUBROUTINE PART_BLOCK(GLOBAL_INDX,N,NP,PV,NV)
IMPLICIT NONE
INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP
INTEGER, INTENT(OUT) :: NV
INTEGER, INTENT(OUT) :: PV(*)
END SUBROUTINE PART_BLOCK
END INTERFACE ! Local variables
INTERFACE
subroutine part_block(global_indx,n,np,pv,nv)
implicit none
integer, intent(in) :: global_indx, n, np
integer, intent(out) :: nv
integer, intent(out) :: pv(*)
end subroutine part_block
end interface ! local variables
interface
! .....user passed subroutine.....
SUBROUTINE PART_BLK2(GLOBAL_INDX,N,NP,PV,NV)
IMPLICIT NONE
INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP
INTEGER, INTENT(OUT) :: NV
INTEGER, INTENT(OUT) :: PV(*)
END SUBROUTINE PART_BLK2
END INTERFACE ! Local variables
subroutine part_blk2(global_indx,n,np,pv,nv)
implicit none
integer, intent(in) :: global_indx, n, np
integer, intent(out) :: nv
integer, intent(out) :: pv(*)
end subroutine part_blk2
end interface ! local variables
INTEGER, PARAMETER :: IZERO=0, IONE=1
CHARACTER, PARAMETER :: ORDER='R'
REAL(KIND(1.D0)), POINTER, SAVE :: B_COL(:), X_COL(:), R_COL(:), &
& B_COL_GLOB(:), X_COL_GLOB(:), R_COL_GLOB(:), B_GLOB(:,:), &
&Z(:), Q(:),Z1(:)
INTEGER :: IARGC, CHECK_DESCR, CONVERT_DESCR
Real(Kind(1.d0)), Parameter :: Dzero = 0.d0, One = 1.d0
Real(Kind(1.d0)) :: MPI_WTIME, T1, T2, TPREC, R_AMAX, B_AMAX,bb(1,1),&
&lambda,scale,resmx,resmxp
integer :: nrhs, nrow, nx1, nx2, n_row, dim,iread
character, parameter :: order='r'
real(kind(1.d0)), pointer, save :: b_col(:), x_col(:), r_col(:), &
& b_col_glob(:), x_col_glob(:), r_col_glob(:)
integer :: iargc
real(kind(1.d0)) :: mpi_wtime, t1, t2, tprec, r_amax, b_amax,bb(1,1),&
&scale,resmx,resmxp
integer :: nrhs, nrow, n_row, dim, nv
logical :: amroot
External IARGC, MPI_WTIME
integer bsze,overlap
common/part/bsze,overlap
INTEGER, POINTER :: WORK(:)
! Sparse Matrices
TYPE(D_SPMAT) :: A, AUX_A, H
TYPE(D_PREC) :: PRE
!!$ TYPE(D_PRECN) :: APRC
! Dense Matrices
REAL(KIND(1.D0)), POINTER :: AUX_B(:,:) , AUX1(:), AUX2(:), VDIAG(:), &
& AUX_G(:,:), AUX_X(:,:), D(:)
integer, pointer :: ivg(:), ipv(:)
! Communications data structure
TYPE(desc_type) :: DESC_A, DESC_A_OUT
external iargc, mpi_wtime
! BLACS parameters
INTEGER :: NPROW, NPCOL, ICTXT, IAM, NP, MYPROW, MYPCOL
! sparse matrices
type(psb_dspmat_type) :: a, aux_a
type(psb_dprec_type) :: pre
! dense matrices
real(kind(1.d0)), pointer :: aux_b(:,:), vdiag(:), d(:)
! communications data structure
type(psb_desc_type):: desc_a, desc_a_out
! blacs parameters
integer :: nprow, npcol, ictxt, iam, np, myprow, mypcol
! solver paramters
integer :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, ml, iprec, novr
! Solver paramters
INTEGER :: ITER, ITMAX, IERR, ITRACE, IRCODE, IPART,&
& METHD, ISTOPC, ML, iprec, novr
integer, pointer :: ierrv(:)
character(len=5) :: afmt
REAL(KIND(1.D0)) :: ERR, EPS
integer iparm(20)
real(kind(1.d0)) rparm(20)
character(len=20) :: name
real(kind(1.d0)) :: err, eps
integer :: iparm(20)
! Other variables
INTEGER :: I,INFO,J
INTEGER :: INTERNAL, M,II,NNZERO
! other variables
integer :: i,info,j
integer :: internal, m,ii,nnzero
! common area
INTEGER M_PROBLEM, NPROC
integer :: m_problem, nproc
allocate(ierrv(6))
! Initialize BLACS
CALL BLACS_PINFO(IAM, NP)
CALL BLACS_GET(IZERO, IZERO, ICTXT)
! initialize blacs
call blacs_pinfo(iam, np)
call blacs_get(izero, izero, ictxt)
! Rectangular Grid, Np x 1
! rectangular grid, np x 1
CALL BLACS_GRIDINIT(ICTXT, ORDER, NP, IONE)
CALL BLACS_GRIDINFO(ICTXT, NPROW, NPCOL, MYPROW, MYPCOL)
AMROOT = (MYPROW==0).AND.(MYPCOL==0)
call blacs_gridinit(ictxt, order, np, ione)
call blacs_gridinfo(ictxt, nprow, npcol, myprow, mypcol)
amroot = (myprow==0).and.(mypcol==0)
name='df_sample'
info=0
!
! Get parameters
! get parameters
!
CALL GET_PARMS(ICTXT,MTRX_FILE,RHS_FILE,CMETHD,PREC,&
& IPART,AFMT,ISTOPC,ITMAX,ITRACE,novr,iprec,EPS)
call get_parms(ictxt,mtrx_file,rhs_file,cmethd,prec,&
& ipart,afmt,istopc,itmax,itrace,novr,iprec,eps)
CALL BLACS_BARRIER(ICTXT,'A')
T1 = MPI_WTIME()
! Read the input matrix to be processed and (possibly) the RHS
NRHS = 1
NPROC = NPROW
call blacs_barrier(ictxt,'a')
t1 = mpi_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
nproc = nprow
IF (AMROOT) THEN
NULLIFY(AUX_B)
CALL READMAT(MTRX_FILE, AUX_A, ICTXT)
if (amroot) then
nullify(aux_b)
call readmat(mtrx_file, aux_a, ictxt)
M_PROBLEM = AUX_A%M
CALL IGEBS2D(ICTXT,'A',' ',1,1,M_PROBLEM,1)
m_problem = aux_a%m
call igebs2d(ictxt,'a',' ',1,1,m_problem,1)
IF(RHS_FILE /= 'NONE') THEN
! Reading an RHS
CALL READ_RHS(RHS_FILE,AUX_B,ICTXT)
END IF
if(rhs_file /= 'NONE') then
! reading an rhs
call read_rhs(rhs_file,aux_b,ictxt)
end if
IF (ASSOCIATED(AUX_B).and.SIZE(AUX_B,1)==M_PROBLEM) THEN
! If any RHS were present, broadcast the first one
write(0,*) 'Ok, got an RHS ',aux_b(m_problem,1)
B_COL_GLOB =>AUX_B(:,1)
ELSE
write(0,*) 'Inventing an RHS '
ALLOCATE(AUX_B(M_PROBLEM,1), STAT=IRCODE)
IF (IRCODE /= 0) THEN
WRITE(0,*) 'Memory allocation failure in DF_SAMPLE'
CALL BLACS_ABORT(ICTXT,-1)
STOP
ENDIF
B_COL_GLOB =>AUX_B(:,1)
DO I=1, M_PROBLEM
!!$ B_COL_GLOB(I) = REAL(I)*2.0/REAL(M_PROBLEM)
B_COL_GLOB(I) = 1.D0
ENDDO
ENDIF
CALL DGEBS2D(ICTXT,'A',' ',M_PROBLEM,1,B_COL_GLOB,M_PROBLEM)
ELSE
CALL IGEBR2D(ICTXT,'A',' ',1,1,M_PROBLEM,1,0,0)
WRITE(0,*) 'Receiving AUX_B'
ALLOCATE(AUX_B(M_PROBLEM,1), STAT=IRCODE)
IF (IRCODE /= 0) THEN
WRITE(0,*) 'Memory allocation failure in DF_SAMPLE'
CALL BLACS_ABORT(ICTXT,-1)
STOP
ENDIF
B_COL_GLOB =>AUX_B(:,1)
CALL DGEBR2D(ICTXT,'A',' ',M_PROBLEM,1,B_COL_GLOB,M_PROBLEM,0,0)
END IF
if (associated(aux_b).and.size(aux_b,1)==m_problem) then
! if any rhs were present, broadcast the first one
write(0,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
write(0,'("Generating an rhs ")')
allocate(aux_b(m_problem,1), stat=ircode)
if (ircode /= 0) then
call psb_errpush(4000,name)
goto 9999
endif
! Switch over different partition types
IF (IPART.EQ.0) THEN
CALL BLACS_BARRIER(ICTXT,'A')
WRITE(6,*) 'Partition type: BLOCK'
CALL MATDIST(AUX_A, A, PART_BLOCK, ICTXT, &
& DESC_A,B_COL_GLOB,B_COL,fmt=afmt)
ELSE IF (IPART.EQ.1) THEN
CALL BLACS_BARRIER(ICTXT,'A')
WRITE(6,*) 'Partition type: BLK2'
CALL MATDIST(AUX_A, A, PART_BLK2, ICTXT, &
& DESC_A,B_COL_GLOB,B_COL,fmt=afmt)
ELSE IF (IPART.EQ.2) THEN
WRITE(0,*) 'Partition type: GRAPH'
IF (AMROOT) THEN
!!$ WRITE(0,*) 'Call BUILD',size(aux_a%ia1),size(aux_a%ia2),np
WRITE(0,*) 'Build type: GRAPH ',aux_a%fida,&
&aux_a%m
CALL BUILD_GRPPART(AUX_A%M,AUX_A%FIDA,AUX_A%IA1,AUX_A%IA2,NP)
ENDIF
CALL BLACS_BARRIER(ICTXT,'A')
!!$ WRITE(0,*) myprow,'Done BUILD_GRPPART'
CALL DISTR_GRPPART(0,0,ICTXT)
!!$ WRITE(0,*) myprow,'Done DISTR_GRPPART'
CALL MATDIST(AUX_A, A, PART_GRAPH, ICTXT, &
& DESC_A,B_COL_GLOB,B_COL)
ELSE
WRITE(6,*) 'Partition type: BLOCK'
CALL MATDIST(AUX_A, A, PART_BLOCK, ICTXT, &
& DESC_A,B_COL_GLOB,B_COL,fmt=afmt)
END IF
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call dgebs2d(ictxt,'a',' ',m_problem,1,b_col_glob,m_problem)
else
call igebr2d(ictxt,'a',' ',1,1,m_problem,1,0,0)
allocate(aux_b(m_problem,1), stat=ircode)
if (ircode /= 0) then
call psb_errpush(4000,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call dgebr2d(ictxt,'a',' ',m_problem,1,b_col_glob,m_problem,0,0)
end if
! switch over different partition types
if (ipart.eq.0) then
call blacs_barrier(ictxt,'a')
write(6,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call matdist(aux_a, a, ivg, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt)
else if (ipart.eq.1) then
call blacs_barrier(ictxt,'a')
write(6,'("Partition type: blk2")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_blk2(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call matdist(aux_a, a, ivg, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt)
else if (ipart.eq.2) then
write(6,'("Partition type: graph")')
if (amroot) then
write(0,'("Build type: graph")')
call build_grppart(aux_a%m,aux_a%fida,aux_a%ia1,aux_a%ia2,np)
endif
call blacs_barrier(ictxt,'a')
call distr_grppart(0,0,ictxt)
call getv_grppart(ivg)
call matdist(aux_a, a, ivg, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt)
else
write(6,'("Partition type: block")')
call matdist(aux_a, a, part_block, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt)
end if
CALL F90_PSDSALL(M_PROBLEM,X_COL,IERRV,DESC_A)
X_COL(:) =0.0
CALL F90_PSDSASB(X_COL,IERRV,DESC_A)
CALL F90_PSDSALL(M_PROBLEM,R_COL,IERRV,DESC_A)
R_COL(:) =0.0
CALL F90_PSDSASB(R_COL,IERRV,DESC_A)
T2 = MPI_WTIME() - T1
call psb_alloc(m_problem,x_col,desc_a,info)
x_col(:) =0.0
call psb_asb(x_col,desc_a,info)
call psb_alloc(m_problem,r_col,desc_a,info)
r_col(:) =0.0
call psb_asb(r_col,desc_a,info)
t2 = mpi_wtime() - t1
CALL DGAMX2D(ICTXT, 'A', ' ', IONE, IONE, T2, IONE,&
& T1, T1, -1, -1, -1)
call dgamx2d(ictxt, 'a', ' ', ione, ione, t2, ione,&
& t1, t1, -1, -1, -1)
IF (AMROOT) THEN
WRITE(6,*) 'Time to Read and Partition Matrix : ',T2
END IF
if (amroot) then
write(6,'("Time to read and partition matrix : ",es10.4)')t2
end if
!
! Prepare the preconditioning matrix. Note the availability
! prepare the preconditioning matrix. note the availability
! of optional parameters
!
IF (AMROOT) WRITE(6,*) 'Preconditioner : "',PREC(1:6),'" ',iprec
if (amroot) write(6,'("Preconditioner : ",a)')prec(1:6)
! Zero initial guess.
! zero initial guess.
select case(iprec)
case(noprec_)
call psb_precset(pre,'noprec')
case(diagsc_)
call psb_precset(pre,'diagsc')
case(ilu_)
case(bja_)
call psb_precset(pre,'ilu')
case(asm_)
call psb_precset(pre,'asm',iv=(/novr,halo_,sum_/))
@@ -217,84 +217,80 @@ PROGRAM DF_SAMPLE
call psb_precset(pre,'asm',iv=(/novr,nohalo_,none_/))
end select
! building the preconditioner
t1 = mpi_wtime()
call psb_precbld(a,pre,desc_a,info)
tprec = mpi_wtime()-t1
call dgamx2d(ictxt,'a',' ',ione, ione,tprec,ione,t1,t1,-1,-1,-1)
write(6,'("Preconditioner time: ",es10.4)')tprec
if (info /= 0) then
write(0,*) 'error in preconditioner :',info
call blacs_abort(ictxt,-1)
stop
end if
iparm = 0
call blacs_barrier(ictxt,'all')
t1 = mpi_wtime()
if (cmethd.eq.'BICGSTAB') then
call psb_bicgstab(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace,istop=istopc)
else if (cmethd.eq.'BICG') then
call psb_bicg(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd.eq.'CGS') then
call psb_cgs(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,itrace)
else if (cmethd.eq.'BICGSTABL') then
call psb_bicgstabl(a,pre,b_col,x_col,eps,desc_a,info,&
& itmax,iter,err,ierr,itrace,ml)
endif
call blacs_barrier(ictxt,'all')
t2 = mpi_wtime() - t1
write(0,*)'Calling gamx2d'
call dgamx2d(ictxt,'a',' ',ione, ione,t2,ione,t1,t1,-1,-1,-1)
write(0,*)'Calling axpby'
call psb_axpby(1.d0,b_col,0.d0,r_col,desc_a,info)
write(0,*)'Calling spmm'
call psb_spmm(-1.d0,a,x_col,1.d0,r_col,desc_a,info)
write(0,*)'Calling nrm2'
call psb_nrm2(resmx,r_col,desc_a,info)
write(0,*)'Calling amax'
call psb_amax(resmxp,r_col,desc_a,info)
T1 = MPI_WTIME()
CALL psb_precbld(A,PRE,DESC_A,INFO)!,'F')
TPREC = MPI_WTIME()-T1
CALL DGAMX2D(ICTXT,'A',' ',IONE, IONE,TPREC,IONE,T1,T1,-1,-1,-1)
WRITE(0,*) 'Preconditioner Time :',TPREC,' ',&
&prec,pre%prec
IF (INFO /= 0) THEN
WRITE(0,*) 'Error in preconditioner :',INFO
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
IPARM = 0
RPARM = 0.D0
CALL BLACS_BARRIER(ICTXT,'All')
T1 = MPI_WTIME()
IF (CMETHD.EQ.'BICGSTAB') Then
CALL F90_BICGSTAB(A,PRE,B_COL,X_COL,EPS,DESC_A,&
& ITMAX,ITER,ERR,IERR,ITRACE,istop=istopc)
!!$ ELSE IF (CMETHD.EQ.'BICG') Then
!!$ CALL F90_BICG(A,IPREC,L,U,VDIAG,B_COL,X_COL,EPS,DESC_A,&
!!$ & ITMAX,ITER,ERR,IERR,ITRACE)
!!$ ELSE IF (CMETHD.EQ.'CGS') Then
!!$ CALL F90_CGS(A,IPREC,L,U,VDIAG,B_COL,X_COL,EPS,DESC_A,&
!!$ & ITMAX,ITER,ERR,IERR,ITRACE)
!!$ ELSE IF (CMETHD.EQ.'BICGSTABL') Then
!!$ CALL F90_BICGSTABL(A,IPREC,L,U,VDIAG,B_COL,X_COL,EPS,DESC_A,&
!!$ & ITMAX,ITER,ERR,IERR,ITRACE,ML)
ENDIF
CALL BLACS_BARRIER(ICTXT,'All')
T2 = MPI_WTIME() - T1
CALL DGAMX2D(ICTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1)
call f90_psaxpby(1.d0,b_col,0.d0,r_col,desc_A)
call f90_psspmm(-1.d0,a,x_col,1.d0,r_col,desc_a)
call f90_nrm2(resmx,r_col,desc_a)
!!$ where (b_col /= 0.d0)
!!$ r_col = r_col/b_col
!!$ end where
call f90_amax(resmxp,r_col,desc_a)
!!$ ITER=IPARM(5)
!!$ ERR = RPARM(2)
!!$ iter=iparm(5)
!!$ err = rparm(2)
if (amroot) then
call prec_descr(6,pre)
call csprt(60+myprow,a)
!!$ write(6,*) 'Number of iterations : ',iter
!!$ WRITE(6,*) 'Error on exit : ',ERR
write(6,*) 'Matrix: ',mtrx_file
write(6,*) 'Computed solution on ',NPROW,' processors.'
write(6,*) 'Iterations to convergence: ',iter
write(6,*) 'Error indicator on exit:',err
write(6,*) 'Time to Buil Prec. : ',TPREC
write(6,*) 'Time to Solve Matrix : ',T2
WRITE(6,*) 'Time per iteration : ',T2/(ITER)
write(6,*) 'Total Time : ',T2+TPREC
write(6,*) 'Residual norm 2 = ',resmx
write(6,*) 'Residual norm inf = ',resmxp
END IF
! call psb_prec_descr(6,pre)
write(6,'("Matrix: ",a)')mtrx_file
write(6,'("Computed solution on ",i4," processors")')nprow
write(6,'("Iterations to convergence: ",i)')iter
write(6,'("Error indicator on exit: ",f7.2)')err
write(6,'("Time to buil prec. : ",es10.4)')tprec
write(6,'("Time to solve matrix : ",es10.4)')t2
write(6,'("Time per iteration : ",es10.4)')t2/(iter)
write(6,'("Total time : ",es10.4)')t2+tprec
write(6,'("Residual norm 2 = ",f7.2)')resmx
write(6,'("Residual norm inf = ",f7.2)')resmxp
end if
allocate(x_col_glob(m_problem),r_col_glob(m_problem),stat=ierr)
if (ierr.ne.0) then
write(0,*) 'Allocation error: no data collection'
write(0,*) 'allocation error: no data collection'
else
call f90_psdgatherm(x_col_glob,x_col,desc_a,iroot=0)
call f90_psdgatherm(r_col_glob,r_col,desc_a,iroot=0)
call psb_gather(x_col_glob,x_col,desc_a,info,iroot=0)
call psb_gather(r_col_glob,r_col,desc_a,info,iroot=0)
if (amroot) then
write(0,*) 'Saving X on file'
write(20,*) 'Matrix: ',mtrx_file
write(20,*) 'Computed solution on ',NPROW,' processors.'
write(20,*) 'Iterations to convergence: ',iter
write(20,*) 'Error indicator (infinity norm) on exit:', &
& ' ||r||/(||A||||x||+||b||) = ',err
write(20,*) 'Max residual = ',resmx, resmxp
write(0,*) 'Saving x on file'
write(20,*) 'matrix: ',mtrx_file
write(20,*) 'computed solution on ',nprow,' processors.'
write(20,*) 'iterations to convergence: ',iter
write(20,*) 'error indicator (infinity norm) on exit:', &
& ' ||r||/(||a||||x||+||b||) = ',err
write(20,*) 'max residual = ',resmx, resmxp
do i=1,m_problem
write(20,998) i,x_col_glob(i),r_col_glob(i),b_col_glob(i)
enddo
@@ -304,15 +300,25 @@ PROGRAM DF_SAMPLE
993 format(i6,4(1x,e12.6))
CALL F90_PSDSFREE(B_COL, DESC_A)
CALL F90_PSDSFREE(X_COL, DESC_A)
CALL F90_PSSPFREE(A, DESC_A)
CALL psb_precfree(PRE,info)
CALL F90_PSDSCFREE(DESC_A,info)
CALL BLACS_GRIDEXIT(ICTXT)
CALL BLACS_EXIT(0)
call psb_free(b_col, desc_a,info)
call psb_free(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_precfree(pre,info)
call psb_dscfree(desc_a,info)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
call blacs_gridexit(ictxt)
call blacs_exit(0)
else
call blacs_gridexit(ictxt)
call blacs_exit(0)
end if
stop
END PROGRAM DF_SAMPLE
end program df_sample
test/Fileread/getp.f90
+1 -1
View File
@@ -123,7 +123,7 @@ CONTAINS
PREC(I:I) = ACHAR(INPARMS(I))
END DO
CALL IGEBR2D(ICONTXT,'A',' ',20,1,INPARMS,20,0,0)
DO I = 1, 20
DO I = 1, LEN(AFMT)
AFMT(I:I) = ACHAR(INPARMS(I))
END DO
test/Fileread/mat_dist.f90
+52 -344
View File
@@ -1,7 +1,7 @@
module mat_dist
interface matdist
module procedure dmatdistf, zmatdistf, dmatdistv
module procedure dmatdistf, dmatdistv
end interface
contains
@@ -64,18 +64,21 @@ contains
! on entry: specifies processor holding a_glob. default: 0
! on exit : unchanged.
!
use f90sparse
implicit none ! parameters
type(d_spmat) :: a_glob
use psb_all_mod
implicit none
! parameters
type(psb_dspmat_type) :: a_glob
real(kind(1.d0)), pointer :: b_glob(:)
integer :: icontxt
type(d_spmat) :: a
type(psb_dspmat_type) :: a
real(kind(1.d0)), pointer :: b(:)
type (desc_type) :: desc_a
type (psb_desc_type) :: desc_a
integer, intent(out) :: info
integer, optional :: inroot
character(len=5), optional :: fmt
interface
! .....user passed subroutine.....
subroutine parts(global_indx,n,np,pv,nv)
implicit none
@@ -83,19 +86,21 @@ contains
integer, intent(out) :: nv
integer, intent(out) :: pv(*)
end subroutine parts
end interface ! local variables
end interface
! local variables
integer :: nprow, npcol, myprow, mypcol
integer :: ircode, length_row, i_count, j_count,&
& k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,&
& i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5)
integer, pointer :: iwork(:)
character :: afmt*5, atyp*5
type(d_spmat) :: blck
type(psb_dspmat_type) :: blck
integer, parameter :: nb=30
real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5, mpi_wtime
external :: mpi_wtime
logical, parameter :: newt=.true.
character(len=20) :: name, ch_err
character(len=20) :: name, ch_err
info = 0
err = 0
@@ -142,7 +147,6 @@ contains
liwork = max(nprow, nrow + ncol)
allocate(iwork(liwork), stat = info)
if (info /= 0) then
write(0,*) 'matdist allocation failed'
info=2025
int_err(1)=liwork
call psb_errpush(info,name,i_err=int_err)
@@ -161,7 +165,7 @@ contains
goto 9999
end if
else
call f90_psdscall(nrow,nrow,parts,icontxt,desc_a,info)
call psb_dscall(nrow,nrow,parts,icontxt,desc_a,info)
if(info/=0) then
info=4010
ch_err='psb_psdscall'
@@ -169,14 +173,14 @@ contains
goto 9999
end if
endif
call f90_psspall(a,desc_a,info,nnz=nnzero/nprow)
call psb_spalloc(a,desc_a,info,nnz=nnzero/nprow)
if(info/=0) then
info=4010
ch_err='psb_psspall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsall(nrow,b,desc_a,info)
call psb_alloc(nrow,b,desc_a,info)
if(info/=0) then
info=4010
ch_err='psb_psdsall'
@@ -188,7 +192,7 @@ contains
blck%m = nb
blck%k = ncol
call spall(blck,nb*ncol,info)
call psb_spall(blck,nb*ncol,info)
if(info/=0) then
info=4010
ch_err='spall'
@@ -201,8 +205,6 @@ contains
do while (i_count.le.nrow)
!!$ write(0,*) myprow,'main loop in matdist',i_count,nrow
!!$ call blacs_barrier(icontxt,'all')
call parts(i_count,nrow,nprow,iwork, length_row)
if (length_row.eq.1) then
@@ -219,8 +221,7 @@ contains
! now we should insert rows i_count..j_count-1
nnr = j_count - i_count
!!$ write(0,*) myprow,'main loop in matdist',i_count,nnr,iproc
!!$ call blacs_barrier(icontxt,'all')
if (myprow == root) then
do j = i_count, j_count
@@ -238,18 +239,18 @@ contains
blck%m = nnr
blck%k = nrow
if (iproc == myprow) then
call f90_psspins(a,i_count,1,blck,desc_a,info)
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
ch_err='psspins'
ch_err='psb_spins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(nnr,b,i_count,b_glob(i_count:j_count-1),&
call psb_ins(nnr,b,i_count,b_glob(i_count:j_count-1),&
&desc_a,info)
if(info/=0) then
info=4010
ch_err='psdsins'
ch_err='psb_ins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
@@ -267,33 +268,33 @@ contains
if (iproc == myprow) then
call igerv2d(icontxt,1,1,nnr,1,root,0)
call igerv2d(icontxt,1,1,ll,1,root,0)
!!$ write(0,*) myprow,'rows and size ',nnr,ll,size(blck%ia2),size(blck%ia1)
call igerv2d(icontxt,nnr+1,1,blck%ia2,nnr+1,root,0)
if (ll > size(blck%ia1)) then
write(0,*) myprow,'need to reallocate ',ll
call spreall(blck,ll,info)
call psb_spreall(blck,ll,info)
if(info/=0) then
info=4010
ch_err='spreall'
ch_err='psb_spreall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
endif
call igerv2d(icontxt,ll,1,blck%ia1,ll,root,0)
call igerv2d(icontxt,nnr+1,1,blck%ia2,nnr+1,root,0)
call dgerv2d(icontxt,ll,1,blck%aspk,ll,root,0)
call dgerv2d(icontxt,nnr,1,b_glob(i_count:i_count+nnr-1),nnr,root,0)
call igesd2d(icontxt,1,1,ll,1,root,0)
blck%m = nnr
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
blck%infoa(psb_nnz_) = ll
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
ch_err='psspins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(nnr,b,i_count,b_glob(i_count:i_count+nnr-1),&
call psb_ins(nnr,b,i_count,b_glob(i_count:i_count+nnr-1),&
&desc_a,info)
if(info/=0) then
info=4010
@@ -330,14 +331,14 @@ contains
blck%m = 1
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
ch_err='psspins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(1,b,i_count,b_glob(i_count:i_count),&
call psb_ins(1,b,i_count,b_glob(i_count:i_count),&
&desc_a,info)
if(info/=0) then
info=4010
@@ -363,14 +364,14 @@ contains
call igesd2d(icontxt,1,1,ll,1,root,0)
blck%m = 1
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
ch_err='psspins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(1,b,i_count,b_glob(i_count:i_count),&
call psb_ins(1,b,i_count,b_glob(i_count:i_count),&
&desc_a,info)
if(info/=0) then
info=4010
@@ -422,7 +423,7 @@ contains
else
call f90_psspasb(a,desc_a,info,dup=1,afmt=afmt)
call psb_spasb(a,desc_a,info,afmt=afmt,dup=1)
if(info/=0)then
info=4010
ch_err='psspasb'
@@ -432,14 +433,14 @@ contains
endif
call f90_psdsasb(b,desc_a,info)
call psb_asb(b,desc_a,info)
if(info/=0)then
info=4010
ch_err='psdsasb'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call spfree(blck,info)
call psb_spfree(blck,info)
if(info/=0)then
info=4010
ch_err='spfree'
@@ -522,14 +523,14 @@ contains
! on entry: specifies processor holding a_glob. default: 0
! on exit : unchanged.
!
use f90sparse
use psb_all_mod
implicit none ! parameters
type(d_spmat) :: a_glob
type(psb_dspmat_type) :: a_glob
real(kind(1.d0)), pointer :: b_glob(:)
integer :: icontxt, v(:)
type(d_spmat) :: a
type(psb_dspmat_type) :: a
real(kind(1.d0)), pointer :: b(:)
type (desc_type) :: desc_a
type (psb_desc_type) :: desc_a
integer, intent(out) :: info
integer, optional :: inroot
character(len=5), optional :: fmt
@@ -540,7 +541,7 @@ contains
& i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5)
integer, pointer :: iwork(:)
character :: afmt*5, atyp*5
type(d_spmat) :: blck
type(psb_dspmat_type) :: blck
integer, parameter :: nb=30
logical, parameter :: newt=.true.
real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5, mpi_wtime
@@ -615,14 +616,14 @@ contains
goto 9999
end if
call f90_psspall(a,desc_a,info,nnz=((nnzero+nprow-1)/nprow))
call psb_spalloc(a,desc_a,info,nnz=((nnzero+nprow-1)/nprow))
if(info/=0) then
info=4010
ch_err='psb_psspall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsall(nrow,b,desc_a,info)
call psb_alloc(nrow,b,desc_a,info)
if(info/=0) then
info=4010
ch_err='psb_psdsall'
@@ -634,7 +635,7 @@ contains
blck%m = nb
blck%k = ncol
call spall(blck,nb*ncol,info)
call psb_spall(blck,nb*ncol,info)
if(info/=0) then
info=4010
ch_err='spall'
@@ -663,7 +664,7 @@ contains
if (myprow == root) then
ll = a_glob%ia2(j_count)-a_glob%ia2(i_count)
if (ll > size(blck%aspk)) then
call spreall(blck,ll,info)
call psb_spreall(blck,ll,info)
if(info/=0) then
info=4010
ch_err='spreall'
@@ -683,17 +684,17 @@ contains
blck%m = nnr
blck%k = nrow
blck%infoa(nnz_) = ll
blck%infoa(psb_nnz_) = ll
if (iproc == myprow) then
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
ch_err='spins'
ch_err='psb_spins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(nnr,b,i_count,b_glob(i_count:j_count-1),&
call psb_ins(nnr,b,i_count,b_glob(i_count:j_count-1),&
&desc_a,info)
if(info/=0) then
info=4010
@@ -717,7 +718,7 @@ contains
call igerv2d(icontxt,1,1,ll,1,root,0)
if (ll > size(blck%aspk)) then
write(0,*) myprow,'need to reallocate ',ll
call spreall(blck,ll,info)
call psb_spreall(blck,ll,info)
if(info/=0) then
info=4010
ch_err='spreall'
@@ -732,7 +733,7 @@ contains
call igesd2d(icontxt,1,1,ll,1,root,0)
blck%m = nnr
blck%k = nrow
blck%infoa(nnz_) = ll
blck%infoa(psb_nnz_) = ll
call psb_spins(ll,blck%ia1,blck%ia2,blck%aspk,a,desc_a,info)
if(info/=0) then
info=4010
@@ -740,7 +741,7 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call f90_psdsins(nnr,b,i_count,b_glob(i_count:i_count+nnr-1),&
call psb_ins(nnr,b,i_count,b_glob(i_count:i_count+nnr-1),&
&desc_a,info)
if(info/=0) then
info=4010
@@ -795,7 +796,7 @@ contains
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call spfree(blck,info)
call psb_spfree(blck,info)
if(info/=0)then
info=4010
ch_err='spfree'
@@ -819,297 +820,4 @@ contains
end subroutine dmatdistv
subroutine zmatdistf (a_glob, a, parts, icontxt, desc_a,&
& b_glob, b, info,inroot,fmt)
!
! an utility subroutine to distribute a matrix among processors
! according to a user defined data distribution, using pessl
! sparse matrix subroutines.
!
! type(d_spmat) :: a_glob
! on entry: this contains the global sparse matrix as follows:
! a%fida =='csr'
! a%aspk for coefficient values
! a%ia1 for column indices
! a%ia2 for row pointers
! a%m for number of global matrix rows
! a%k for number of global matrix columns
! on exit : undefined, with unassociated pointers.
!
! type(d_spmat) :: a
! on entry: fresh variable.
! 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, intent(in) :: global_indx, n, np
! integer, intent(out) :: nv
! integer, 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 :: icontxt
! on entry: blacs context.
! on exit : unchanged.
!
! type (desc_type) :: desc_a
! on entry: fresh variable.
! on exit : the updated array descriptor
!
! real(kind(1.d0)), pointer, optional :: b_glob(:)
! on entry: this contains right hand side.
! on exit :
!
! real(kind(1.d0)), pointer, optional :: b(:)
! on entry: fresh variable.
! on exit : this will contain the local right hand side.
!
! integer, optional :: inroot
! on entry: specifies processor holding a_glob. default: 0
! on exit : unchanged.
!
use typesp
use typedesc
use f90tools
implicit none ! parameters
type(z_spmat) :: a_glob
complex(kind(1.d0)), pointer :: b_glob(:)
integer :: icontxt
type(z_spmat) :: a
complex(kind(1.d0)), pointer :: b(:)
type (desc_type) :: desc_a
integer, intent(out) :: info
integer, optional :: inroot
character(len=5), optional :: fmt
interface
! .....user passed subroutine.....
subroutine parts(global_indx,n,np,pv,nv)
implicit none
integer, intent(in) :: global_indx, n, np
integer, intent(out) :: nv
integer, intent(out) :: pv(*)
end subroutine parts
end interface ! local variables
integer :: nprow, npcol, myprow, mypcol
integer :: ircode, length_row, i_count, j_count,&
& k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,&
& i,j,k, ll, isize, iproc, nnr, err, err_act, int_err(5)
integer, pointer :: iwork(:)
character :: afmt*5, atyp*5
type(z_spmat) :: blck
integer, parameter :: nb = 30
character(len=20) :: name, ch_err
info = 0
err = 0
name = 'mat_distf'
call psb_erractionsave(err_act)
! executable statements
if (present(inroot)) then
root = inroot
else
root = 0
end if
call blacs_gridinfo(icontxt, nprow, npcol, myprow, mypcol)
if (myprow == root) then
! extract information from a_glob
if (a_glob%fida.ne. 'CSR') then
write(0,*) 'unsupported input matrix format'
call blacs_abort(icontxt,-1)
endif
nrow = a_glob%m
ncol = a_glob%k
if (nrow /= ncol) then
write(0,*) 'a rectangular matrix ? ',nrow,ncol
call blacs_abort(icontxt,-1)
endif
nnzero = size(a_glob%aspk)
nrhs = 1
! broadcast informations to other processors
call igebs2d(icontxt, 'a', ' ', 1, 1, nrow, 1)
call igebs2d(icontxt, 'a', ' ', 1, 1, ncol, 1)
call igebs2d(icontxt, 'a', ' ', 1, 1, nnzero, 1)
call igebs2d(icontxt, 'a', ' ', 1, 1, nrhs, 1)
else !(myprow /= root)
! receive informations
call igebr2d(icontxt, 'a', ' ', 1, 1, nrow, 1, root, 0)
call igebr2d(icontxt, 'a', ' ', 1, 1, ncol, 1, root, 0)
call igebr2d(icontxt, 'a', ' ', 1, 1, nnzero, 1, root, 0)
call igebr2d(icontxt, 'a', ' ', 1, 1, nrhs, 1, root, 0)
end if ! allocate integer work area
liwork = max(nprow, nrow + ncol)
allocate(iwork(liwork), stat = ircode)
if (ircode /= 0) then
write(0,*) 'matdist allocation failed'
return
endif
if (myprow == root) then
write (*, fmt = *) 'start matdist',root, size(iwork)
endif
call f90_psdscall(nrow,nrow,parts,icontxt,desc_a,info)
call f90_psspall(a,desc_a,info,nnz=nnzero/nprow)
call f90_psdsall(nrow,b,desc_a,info)
isize = max(3*nb,ncol)
allocate(blck%aspk(nnzero),blck%ia1(nnzero),blck%ia2(nnzero),stat=ircode)
if (ircode /= 0) then
write(0,*) 'error on allocating blck'
call blacs_abort(icontxt,-1)
stop
endif
blck%m = 1
blck%k = ncol
blck%fida = 'csr'
i_count = 1
do while (i_count.le.nrow)
call parts(i_count,nrow,nprow,iwork, length_row)
if (length_row.eq.1) then
j_count = i_count + 1
iproc = iwork(1)
call parts(j_count,nrow,nprow,iwork, length_row)
do while ((j_count.le.nrow).and.(j_count-i_count.lt.nb)&
&.and.(length_row.eq.1).and.(iwork(1).eq.iproc))
j_count = j_count + 1
if (j_count.le.nrow) &
& call parts(j_count,nrow,nprow,iwork, length_row)
end do
! now we should insert rows i_count..j_count-1
nnr = j_count - i_count
if (myprow == root) then
do j = i_count, j_count
blck%ia2(j-i_count+1) = a_glob%ia2(j) - &
& a_glob%ia2(i_count) + 1
enddo
k = a_glob%ia2(i_count)
do j = k, a_glob%ia2(j_count)-1
blck%aspk(j-k+1) = a_glob%aspk(j)
blck%ia1(j-k+1) = a_glob%ia1(j)
enddo
ll = blck%ia2(nnr+1) - 1
blck%m = nnr
blck%k = nrow
if (iproc == myprow) then
call f90_psspins(a,i_count,1,blck,desc_a,info)
call f90_psdsins(nnr,b,i_count,b_glob(i_count:j_count-1),&
&desc_a,info)
else
call igesd2d(icontxt,1,1,nnr,1,iproc,0)
call igesd2d(icontxt,1,1,ll,1,iproc,0)
call igesd2d(icontxt,nnr+1,1,blck%ia2,nnr+1,iproc,0)
call igesd2d(icontxt,ll,1,blck%ia1,ll,iproc,0)
call zgesd2d(icontxt,ll,1,blck%aspk,ll,iproc,0)
call zgesd2d(icontxt,nnr,1,b_glob(i_count:j_count-1),nnr,iproc,0)
call igerv2d(icontxt,1,1,ll,1,iproc,0)
endif
else if (myprow /= root) then
if (iproc == myprow) then
call igerv2d(icontxt,1,1,nnr,1,root,0)
call igerv2d(icontxt,1,1,ll,1,root,0)
call igerv2d(icontxt,nnr+1,1,blck%ia2,nnr+1,root,0)
call igerv2d(icontxt,ll,1,blck%ia1,ll,root,0)
call zgerv2d(icontxt,ll,1,blck%aspk,ll,root,0)
call zgerv2d(icontxt,nnr,1,b_glob(i_count:i_count+nnr-1),nnr,root,0)
call igesd2d(icontxt,1,1,ll,1,root,0)
blck%m = nnr
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
call f90_psdsins(nnr,b,i_count,b_glob(i_count:i_count+nnr-1),&
&desc_a,info)
endif
endif
i_count = j_count
else
! here processors are counted 1..nprow
do j_count = 1, length_row
k_count = iwork(j_count)
if (myprow == root) then
blck%ia2(1) = 1
blck%ia2(2) = 1
do j = a_glob%ia2(i_count), a_glob%ia2(i_count+1)-1
blck%aspk(blck%ia2(2)) = a_glob%aspk(j)
blck%ia1(blck%ia2(2)) = a_glob%ia1(j)
blck%ia2(2) =blck%ia2(2) + 1
enddo
ll = blck%ia2(2) - 1
if (k_count == myprow) then
blck%infoa(1) = ll
blck%infoa(2) = ll
blck%infoa(3) = 2
blck%infoa(4) = 1
blck%infoa(5) = 1
blck%infoa(6) = 1
blck%m = 1
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
call f90_psdsins(1,b,i_count,b_glob(i_count:i_count),&
&desc_a,info)
else
call igesd2d(icontxt,1,1,ll,1,k_count,0)
call igesd2d(icontxt,ll,1,blck%ia1,ll,k_count,0)
call zgesd2d(icontxt,ll,1,blck%aspk,ll,k_count,0)
call zgesd2d(icontxt,1,1,b_glob(i_count),1,k_count,0)
call igerv2d(icontxt,1,1,ll,1,k_count,0)
endif
else if (myprow /= root) then
if (k_count == myprow) then
call igerv2d(icontxt,1,1,ll,1,root,0)
blck%ia2(1) = 1
blck%ia2(2) = ll+1
call igerv2d(icontxt,ll,1,blck%ia1,ll,root,0)
call zgerv2d(icontxt,ll,1,blck%aspk,ll,root,0)
call zgerv2d(icontxt,1,1,b_glob(i_count),1,root,0)
call igesd2d(icontxt,1,1,ll,1,root,0)
blck%m = 1
blck%k = nrow
call f90_psspins(a,i_count,1,blck,desc_a,info)
call f90_psdsins(1,b,i_count,b_glob(i_count:i_count),&
&desc_a,info)
endif
endif
end do
i_count = i_count + 1
endif
end do
! default storage format for sparse matrix; we do not
! expect duplicated entries.
if (present(fmt)) then
afmt=fmt
else
afmt = 'csr'
endif
call f90_psspasb(a,desc_a,info,dup=1)
call f90_psdsasb(b,desc_a,info)
call spfree(blck,info)
deallocate(iwork)
if (myprow == root) write (*, fmt = *) 'end matdist'
return
end subroutine zmatdistf
end module mat_dist
test/Fileread/mmio.f90
+9 -259
View File
@@ -1,20 +1,20 @@
module mmio
use typesp
use psb_spmat_type
public mm_mat_read, mm_mat_write
interface mm_mat_read
module procedure dmm_mat_read, zmm_mat_read
module procedure dmm_mat_read
end interface
interface mm_mat_write
module procedure dmm_mat_write
end interface
private desym,zdesym
private desym
contains
subroutine dmm_mat_read(a, iret, iunit, filename)
use typesp
use psb_spmat_type
implicit none
type(d_spmat), intent(out) :: a
type(psb_dspmat_type), intent(out) :: a
integer, intent(out) :: iret
integer, optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename
@@ -140,7 +140,7 @@ contains
call desym(nrow, a%aspk, a%ia2, a%ia1, as_loc, ia2_loc,&
& ia1_loc, iwork, nnzero, nzr)
call spreall(a,nzr,ircode)
call psb_spreall(a,nzr,ircode)
if (ircode /= 0) goto 993
allocate(tmp(nzr),stat=ircode)
if (ircode /= 0) goto 993
@@ -189,178 +189,13 @@ contains
end subroutine dmm_mat_read
subroutine zmm_mat_read(a, iret, iunit, filename)
use typesp
implicit none
type(z_spmat), intent(out) :: a
integer, intent(out) :: iret
integer, optional, intent(in) :: iunit
character(len=*), optional, intent(in) :: filename
character :: mmheader*15, fmt*15, object*10, type*10, sym*15, line*1024
integer :: indcrd, ptrcrd, totcrd,&
& valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
complex(kind(1.0d0)), pointer :: as_loc(:), dwork(:)
integer, pointer :: ia1_loc(:), ia2_loc(:), iwork(:), tmp(:), aux(:)
integer :: ircode, i,iel,ptr,nzr,infile,&
& j, liwork, ldwork, root, nprow, npcol, myprow, mypcol
iret = 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=*,end=902) mmheader, object, fmt, type, sym
call lowerc(object,1,10)
call lowerc(fmt,1,15)
if ( (object .ne. 'matrix').or.(fmt.ne.'coordinate')) then
write(0,*) 'READ_MATRIX: input file type not yet supported'
iret=909
return
end if
do
read(infile,fmt='(a)') line
if (line(1:1) /= '%') exit
end do
read(line,fmt=*) nrow,ncol,nnzero
a%m = nrow
a%k = ncol
a%fida = 'CSR'
call lowerc(type,1,10)
call lowerc(sym,1,15)
if ((type == 'complex').and.(sym == 'general')) then
allocate(a%aspk(nnzero), a%ia1(nnzero), a%ia2(nrow+1),&
& a%pl(nrow),a%pr(nrow), tmp(nnzero+1), aux(nnzero+2),stat = ircode)
if (ircode /= 0) goto 993
do i=1,nnzero
read(infile,fmt=*,end=902) tmp(i),a%ia1(i),a%aspk(i)
end do
call mrgsrt(nnzero,tmp,aux,ircode)
if (ircode.eq.0) call zreordvn(nnzero,a%aspk,tmp,a%ia1,aux)
! .... Order with key a%ia1 (COLUMN INDEX) ...
i = 1
j = i
! .... order with key tmp (row index) ...
do
if (i > nnzero) exit
do
if (j > nnzero) exit
if (tmp(j) /= tmp(i)) exit
j = j+1
! if (j.eq.(nnzero+1)) exit
enddo
iel = j - i
call mrgsrt(iel,a%ia1(i),aux,ircode)
if (ircode == 0) call zreordvn(iel,a%aspk(i),tmp(i),&
& a%ia1(i), aux)
i = j
enddo
! convert to csr format
iel = 1
a%ia2(1) = 1
do i = a%ia2(1), nrow
do
if (iel > nnzero) exit
if (tmp(iel) /= i) exit
iel = iel + 1
enddo
a%ia2(i+1) = iel
enddo
deallocate(aux,tmp)
else if ((type == 'complex').and.(sym == 'symmetric')) then
! we are generally working with non-symmetric matrices, so
! we de-symmetrize what we are about to read
allocate(a%aspk(2*nnzero),a%ia1(2*nnzero),&
& a%ia2(2*nnzero),as_loc(2*nnzero),&
& ia1_loc(2*nnzero),ia2_loc(2*nnzero),&
&a%pl(nrow),a%pr(nrow), stat = ircode)
if (ircode /= 0) goto 993
do i=1,nnzero
read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),a%aspk(i)
end do
liwork = 2*nnzero+2
allocate(iwork(liwork), stat = ircode)
if (ircode /= 0) goto 993
! After this call NNZERO contains the actual value for
! desymetrized matrix
call zdesym(nrow, a%aspk, a%ia2, a%ia1, as_loc, ia2_loc,&
& ia1_loc, iwork, nnzero, nzr)
deallocate(a%aspk,a%ia1,a%ia2)
nnzero=nzr
!!$ call spreall(a,nzr,ircode)
if (ircode /= 0) goto 993
allocate(tmp(nzr),stat=ircode)
if (ircode /= 0) goto 993
a%aspk(1:nzr) = as_loc(1:nzr)
a%ia1(1:nzr) = ia2_loc(1:nzr)
tmp(1:nzr) = ia1_loc(1:nzr)
iel = 1
a%ia2(1) = 1
do i = 1, nrow
do
if (tmp(iel) /= i) exit
iel = iel + 1
if (iel > nzr) exit
enddo
a%ia2(i+1) = iel
enddo
deallocate(as_loc, ia1_loc, ia2_loc,tmp,iwork)
else
write(0,*) 'read_matrix: matrix type not yet supported'
iret=904
end if
if (infile/=5) close(infile)
return
! open failed
901 iret=901
write(0,*) 'read_matrix: could not open file ',filename,' for input'
return
902 iret=902
write(0,*) 'READ_MATRIX: Unexpected end of file '
return
993 iret=993
write(0,*) 'READ_MATRIX: Memory allocation failure'
return
end subroutine zmm_mat_read
subroutine dmm_mat_write(a,mtitle,iret,eiout,filename)
use typesp
use psb_spmat_type
implicit none
type(d_spmat), intent(in) :: a
type(psb_dspmat_type), intent(in) :: a
integer, intent(out) :: iret
character(len=*), intent(in) :: mtitle
integer, optional, intent(in) :: eiout
@@ -389,28 +224,11 @@ contains
endif
endif
call dcsprt(a%m,a%k,a%fida,a%descra,a%aspk,a%ia1,a%ia2,a%infoa,&
& mtitle,iout,iret)
call psb_dcsprt(iout,a)
if (iout /= 6) close(iout)
!!$ write(outfile(9:),998) '.xrhs'
!!$ open (iout,file=outfile,status='replace',err=901)
!!$ write(iout,fmt=997)
!!$ write(iout,fmt=996) mtitle
!!$ write(iout,fmt=995) 'Number of equations ',nrow
!!$ write(iout,fmt=995) 'Number of iterations to convergence ',iter
!!$ write(iout,fmt=996)
!!$ write(iout,fmt=996) 'index comp. solution Right hand side'
!!$ write(iout,fmt=997)
!!$ do i=1, nrow
!!$ write(iout,993) i,x(i),rhs(i)
!!$993 format(i5,4(1x,e12.6))
!!$ enddo
!!$ close(iout)
!!$ !$$$ call system('gzip -f9 '//outfile)
return
901 continue
@@ -420,20 +238,6 @@ contains
end subroutine dmm_mat_write
!!$ subroutine lowerc(string,pos,len)
!!$ integer pos, len
!!$ character(len=*) string
!!$
!!$ character(len=26), parameter :: lcase='abcdefghijklmnopqrstuvwxyz',&
!!$ & ucase='ABCDEFGHIJKLMNOPQRSTUVWXYZ'
!!$
!!$ do i=pos,pos+len-1
!!$ k = index(ucase,string(i:i))
!!$ if (k.ne.0) string(i:i) = lcase(k:k)
!!$ enddo
!!$ return
!!$ end subroutine lowerc
subroutine desym(nrow,a,ja,ia,as,jas,ias,aux,nnzero,nzr)
implicit none
! .. scalar arguments ..
@@ -491,58 +295,4 @@ contains
return
end subroutine desym
subroutine zdesym(nrow,a,ja,ia,as,jas,ias,aux,nnzero,nzr)
implicit none
! .. scalar arguments ..
integer :: nrow,nnzero,value,index,ptr, nzr
! .. array arguments ..
complex(kind(1.d0)) :: a(*),as(*)
integer :: ia(*),ias(*),jas(*),ja(*),aux(*)
! .. local scalars ..
integer :: i,iaw1,iaw2,iawt,j,jpt,k,kpt,ldim,nzl,js,iret,nel,diagel
! ..
nel = 0
diagel=0
do i=1, nnzero
as(i) = a(i)
jas(i) = ja(i)
ias(i) = ia(i)
if(ja(i) < ia(i)) then !this control avoids malfunctions in the cases
! where the matrix is declared symmetric but all its elements are
! explicitly stored see young1c.mtx from "Matrix Market".
! Nominally Matrix Market only stores lower triangle.
nel = nel+1
as(nnzero+nel) = a(i)
jas(nnzero+nel) = ia(i)
ias(nnzero+nel) = ja(i)
end if
end do
! .... order with key ias ...
nzr = nnzero + nel
call mrgsrt(nzr,ias,aux,iret)
if (iret == 0) call zreordvn(nzr,as,ias,jas,aux)
! .... order with key jas ...
i = 1
j = i
do
if (i > nzr) exit
do
if (j > nzr) exit
if (ias(j) /= ias(i)) exit
j = j+1
enddo
nzl = j - i
call mrgsrt(nzl,jas(i),aux,iret)
if (iret.eq.0) call zreordvn(nzl,as(i),ias(i),jas(i),aux)
i = j
enddo
return
end subroutine zdesym
end module mmio
test/Fileread/partgraph.f90
-1
View File
@@ -104,7 +104,6 @@ CONTAINS
SUBROUTINE BUILD_GRPPART(N,FIDA,IA1,IA2,NPARTS)
USE TYPESP
INTEGER :: NPARTS
INTEGER :: IA1(:), IA2(:)
INTEGER :: N, I, IB, II,numflag,nedc,wgflag
test/Fileread/read_mat.f90
+50 -138
View File
@@ -43,11 +43,11 @@ module read_mat
contains
subroutine readmat (filename, a, ictxt, inroot)
use typesp
use psb_spmat_type
use mmio
implicit none
integer :: ictxt
type(d_spmat) :: a
type(psb_dspmat_type) :: a
character(len=*) :: filename
integer, optional :: inroot
integer, parameter :: infile = 2
@@ -72,17 +72,20 @@ contains
end subroutine readmat
subroutine zreadmat (filename, a, ictxt, inroot)
use typesp
use mmio
subroutine read_rhs (filename, b, ictxt, inroot,&
& indwork, iniwork)
implicit none
integer :: ictxt
type(z_spmat) :: a
character :: filename*(*)
integer, optional :: inroot
integer, parameter :: infile = 2
integer :: info, root, nprow, npcol, myprow, mypcol
real(kind(1.0d0)), pointer, optional :: indwork(:)
integer, pointer, optional :: iniwork(:) ! local variables
integer, parameter :: infile = 2
integer :: nrow, ncol, i,root, nprow, npcol, myprow, mypcol, ircode, j
character :: mmheader*15, fmt*15, object*10, type*10, sym*15,&
& line*1024
real(kind(1.0d0)), pointer :: b(:,:)
if (present(inroot)) then
root = inroot
else
@@ -90,138 +93,47 @@ contains
end if
call blacs_gridinfo(ictxt, nprow, npcol, myprow, mypcol)
if (myprow == root) then
write(*, *) 'start read_matrix' ! open input file
call mm_mat_read(a,info,infile,filename)
if (info /= 0) then
write(0,*) 'Error return from MM_MAT_READ ',info
call blacs_abort(ictxt, 1) ! Unexpected End of File
endif
write(*, *) 'start read_rhs' ! open input file
open(infile,file=filename, status='old', err=901, action="read")
read(infile,fmt=*, end=902) mmheader, object, fmt, type, sym
write(0,*)'obj fmt',object, fmt
if ( (object .ne. 'matrix').or.(fmt.ne.'array')) then
write(0,*) 'read_rhs: input file type not yet supported'
call blacs_abort(ictxt, 1)
end if
do
read(infile,fmt='(a)') line
if (line(1:1) /= '%') exit
end do
read(line,fmt=*)nrow,ncol
call lowerc(type,1,10)
call lowerc(sym,1,15)
if ((type == 'real').and.(sym == 'general')) then
allocate(b(nrow,ncol),stat = ircode)
if (ircode /= 0) goto 993
read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol)
else
write(0,*) 'read_rhs: rhs type not yet supported'
call blacs_abort(ictxt, 1)
end if ! read right hand sides
write(*,*) 'end read_rhs'
end if
return
end subroutine zreadmat
SUBROUTINE READ_RHS (FILENAME, B, ICTXT, INROOT,&
& INDWORK, INIWORK)
IMPLICIT NONE
INTEGER :: ICTXT
CHARACTER :: FILENAME*(*)
INTEGER, OPTIONAL :: INROOT
REAL(KIND(1.0D0)), POINTER, OPTIONAL :: INDWORK(:)
INTEGER, POINTER, OPTIONAL :: INIWORK(:) ! Local Variables
INTEGER, PARAMETER :: INFILE = 2
INTEGER :: NROW, NCOL, I,ROOT, NPROW, NPCOL, MYPROW, MYPCOL, IRCODE, J
CHARACTER :: MMHEADER*15, FMT*15, OBJECT*10, TYPE*10, SYM*15,&
& LINE*1024
REAL(KIND(1.0D0)), POINTER :: B(:,:)
IF (PRESENT(INROOT)) THEN
ROOT = INROOT
ELSE
ROOT = 0
END IF
CALL BLACS_GRIDINFO(ICTXT, NPROW, NPCOL, MYPROW, MYPCOL)
IF (MYPROW == ROOT) THEN
WRITE(*, *) 'Start read_rhs' ! Open Input File
OPEN(INFILE,FILE=FILENAME, STATUS='OLD', ERR=901, ACTION="READ")
READ(INFILE,FMT=*, END=902) MMHEADER, OBJECT, FMT, TYPE, SYM
write(0,*)'obj fmt',object, fmt
IF ( (OBJECT .NE. 'matrix').OR.(FMT.NE.'array')) THEN
WRITE(0,*) 'READ_RHS: input file type not yet supported'
CALL BLACS_ABORT(ICTXT, 1)
END IF
do
read(infile,fmt='(a)') line
if (line(1:1) /= '%') exit
end do
READ(LINE,FMT=*)NROW,NCOL
CALL LOWERC(TYPE,1,10)
CALL LOWERC(SYM,1,15)
IF ((TYPE == 'real').AND.(SYM == 'general')) THEN
ALLOCATE(B(NROW,NCOL),STAT = IRCODE)
IF (IRCODE /= 0) GOTO 993
READ(INFILE,FMT=*,END=902) ((B(I,J), I=1,NROW),J=1,NCOL)
ELSE
WRITE(0,*) 'READ_RHS: RHS type not yet supported'
CALL BLACS_ABORT(ICTXT, 1)
END IF ! Read Right Hand Sides
WRITE(*,*) 'End READ_RHS'
END IF
RETURN
! Open failed
901 WRITE(0,*) 'READ_RHS: Could not open file ',&
& INFILE,' for input'
CALL BLACS_ABORT(ICTXT, 1) ! Unexpected End of File
902 WRITE(0,*) 'READ_RHS: Unexpected end of file ',INFILE,&
! open failed
901 write(0,*) 'read_rhs: could not open file ',&
& infile,' for input'
call blacs_abort(ictxt, 1) ! unexpected end of file
902 write(0,*) 'read_rhs: unexpected end of file ',infile,&
& ' during input'
CALL BLACS_ABORT(ICTXT, 1) ! Allocation Failed
993 WRITE(0,*) 'READ_RHS: Memory allocation failure'
CALL BLACS_ABORT(ICTXT, 1)
END SUBROUTINE READ_RHS
SUBROUTINE ZREAD_RHS(FILENAME, B, ICTXT, INROOT)
IMPLICIT NONE
INTEGER :: ICTXT
CHARACTER :: FILENAME*(*)
INTEGER, OPTIONAL :: INROOT
INTEGER, PARAMETER :: INFILE = 2
INTEGER :: NROW, NCOL, I,ROOT, NPROW, NPCOL, MYPROW, MYPCOL, IRCODE, J
CHARACTER :: MMHEADER*15, FMT*15, OBJECT*10, TYPE*10, SYM*15,&
& LINE*1024
COMPLEX(KIND(1.0D0)), POINTER :: B(:,:)
IF (PRESENT(INROOT)) THEN
ROOT = INROOT
ELSE
ROOT = 0
END IF
CALL BLACS_GRIDINFO(ICTXT, NPROW, NPCOL, MYPROW, MYPCOL)
IF (MYPROW == ROOT) THEN
WRITE(*, *) 'Start read_rhs' ! Open Input File
OPEN(INFILE,FILE=FILENAME, STATUS='OLD', ERR=901, ACTION="READ")
READ(INFILE,FMT=*, END=902) MMHEADER, OBJECT, FMT, TYPE, SYM
write(0,*)'obj fmt',object, fmt
IF ( (OBJECT .NE. 'matrix').OR.(FMT.NE.'array')) THEN
WRITE(0,*) 'READ_RHS: input file type not yet supported'
CALL BLACS_ABORT(ICTXT, 1)
END IF
do
read(infile,fmt='(a)') line
if (line(1:1) /= '%') exit
end do
READ(LINE,FMT=*)NROW,NCOL
CALL LOWERC(TYPE,1,10)
CALL LOWERC(SYM,1,15)
IF ((TYPE == 'complex').AND.(SYM == 'general')) THEN
ALLOCATE(B(NROW,NCOL),STAT = IRCODE)
IF (IRCODE /= 0) GOTO 993
READ(INFILE,FMT=*,END=902) ((B(I,J), I=1,NROW),J=1,NCOL)
ELSE
WRITE(0,*) 'READ_RHS: RHS type not yet supported'
CALL BLACS_ABORT(ICTXT, 1)
END IF ! Read Right Hand Sides
WRITE(*,*) 'End READ_RHS'
END IF
RETURN
! Open failed
901 WRITE(0,*) 'READ_RHS: Could not open file ',&
& INFILE,' for input'
CALL BLACS_ABORT(ICTXT, 1) ! Unexpected End of File
902 WRITE(0,*) 'READ_RHS: Unexpected end of file ',INFILE,&
& ' during input'
CALL BLACS_ABORT(ICTXT, 1) ! Allocation Failed
993 WRITE(0,*) 'READ_RHS: Memory allocation failure'
CALL BLACS_ABORT(ICTXT, 1)
END SUBROUTINE ZREAD_RHS
call blacs_abort(ictxt, 1) ! allocation failed
993 write(0,*) 'read_rhs: memory allocation failure'
call blacs_abort(ictxt, 1)
end subroutine read_rhs
END MODULE READ_MAT
end module read_mat