From 0383e806183c7985ee385754bcf5c3a4d036ae0b Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 12 Jan 2007 09:13:36 +0000 Subject: [PATCH] Reorganization of test dir. --- test/Fileread/Makefile | 66 - test/Fileread/hbio.f90 | 638 ------- test/Fileread/mat_dist.f90 | 1584 ----------------- test/Fileread/mmio.f90 | 379 ---- test/Fileread/part_blk2.f | 86 - test/Fileread/part_block.f | 97 - test/Fileread/partgraph.f90 | 222 --- test/Fileread/read_mat.f90 | 253 --- test/fileread/Makefile | 36 + test/{Fileread => fileread}/df_sample.f90 | 53 +- test/{Fileread => fileread}/getp.f90 | 2 +- .../{Fileread/RUNS => fileread/runs}/Makefile | 0 test/{Fileread/RUNS => fileread/runs}/dfs.inp | 0 test/{Fileread/RUNS => fileread/runs}/zfs.inp | 0 test/{Fileread => fileread}/zf_sample.f90 | 54 +- test/pargen/Makefile | 14 +- test/pargen/part_block.f | 80 - test/pargen/ppde90.f90 | 23 +- test/pargen/{RUNS => runs}/Makefile | 0 test/pargen/{RUNS => runs}/mach | 0 test/pargen/{RUNS => runs}/ppde.inp | 2 +- test/util/Makefile | 42 + test/{Fileread => util}/dhb2mm.f90 | 5 +- test/{Fileread => util}/dmm2hb.f90 | 5 +- test/{Fileread => util}/zhb2mm.f90 | 5 +- test/{Fileread => util}/zmm2hb.f90 | 5 +- 26 files changed, 114 insertions(+), 3537 deletions(-) delete mode 100644 test/Fileread/Makefile delete mode 100644 test/Fileread/hbio.f90 delete mode 100644 test/Fileread/mat_dist.f90 delete mode 100644 test/Fileread/mmio.f90 delete mode 100644 test/Fileread/part_blk2.f delete mode 100644 test/Fileread/part_block.f delete mode 100644 test/Fileread/partgraph.f90 delete mode 100644 test/Fileread/read_mat.f90 create mode 100644 test/fileread/Makefile rename test/{Fileread => fileread}/df_sample.f90 (85%) rename test/{Fileread => fileread}/getp.f90 (99%) rename test/{Fileread/RUNS => fileread/runs}/Makefile (100%) rename test/{Fileread/RUNS => fileread/runs}/dfs.inp (100%) rename test/{Fileread/RUNS => fileread/runs}/zfs.inp (100%) rename test/{Fileread => fileread}/zf_sample.f90 (85%) delete mode 100644 test/pargen/part_block.f rename test/pargen/{RUNS => runs}/Makefile (100%) rename test/pargen/{RUNS => runs}/mach (100%) rename test/pargen/{RUNS => runs}/ppde.inp (87%) create mode 100644 test/util/Makefile rename test/{Fileread => util}/dhb2mm.f90 (98%) rename test/{Fileread => util}/dmm2hb.f90 (97%) rename test/{Fileread => util}/zhb2mm.f90 (98%) rename test/{Fileread => util}/zmm2hb.f90 (97%) diff --git a/test/Fileread/Makefile b/test/Fileread/Makefile deleted file mode 100644 index 5a8c4150..00000000 --- a/test/Fileread/Makefile +++ /dev/null @@ -1,66 +0,0 @@ -include ../../Make.inc -# -# Libraries used -# -LIBDIR=../../lib/ -PSBLAS_LIB= -L$(LIBDIR) -lpsblas - -# -# We are using the public domain tool METIS from U. Minnesota. To get it -# check URL http://www.cs.umn.edu:~karypis -# -METIS_LIB = -L$(HOME)/NUMERICAL/metis-4.0 -lmetis - -INCDIRS=-I$(LIBDIR) - -DFOBJS=partgraph.o part_block.o read_mat.o getp.o \ - mmio.o mat_dist.o df_sample.o part_blk2.o -ZFOBJS=partgraph.o part_block.o read_mat.o getp.o \ - mmio.o mat_dist.o zf_sample.o part_blk2.o -IOOBJS= mmio.o hbio.o - -ZH2MOBJS=zhb2mm.o $(IOOBJS) -DH2MOBJS=dhb2mm.o $(IOOBJS) -DM2HOBJS=dmm2hb.o $(IOOBJS) -ZM2HOBJS=zmm2hb.o $(IOOBJS) -MMHBOBJS=zhb2mm.o dhb2mm.o dmm2hb.o zmm2hb.o - -EXEDIR=./RUNS - -all: df_sample zf_sample dhb2mm zhb2mm dmm2hb zmm2hb - -read_mat.o: mmio.o - -df_sample: $(DFOBJS) - $(F90LINK) $(LINKOPT) $(DFOBJS) -o df_sample\ - $(PSBLAS_LIB) $(METIS_LIB) $(BLACS) $(SLU) $(UMF) $(BLAS) - /bin/mv df_sample $(EXEDIR) -zf_sample: $(ZFOBJS) - $(F90LINK) $(LINKOPT) $(ZFOBJS) -o zf_sample\ - $(PSBLAS_LIB) $(METIS_LIB) $(BLACS) $(SLU) $(UMF) $(BLAS) - /bin/mv zf_sample $(EXEDIR) - -$(MMHBOBJS): $(IOOBJS) -dhb2mm: $(DH2MOBJS) - $(MPF90) -o dhb2mm $(DH2MOBJS) $(PSBLAS_LIB) $(BLACS) -dmm2hb: $(DM2HOBJS) - $(MPF90) -o dmm2hb $(DM2HOBJS) $(PSBLAS_LIB) $(BLACS) -zhb2mm: $(ZH2MOBJS) - $(MPF90) -o zhb2mm $(ZH2MOBJS) $(PSBLAS_LIB) $(BLACS) -zmm2hb: $(ZM2HOBJS) - $(MPF90) -o zmm2hb $(ZM2HOBJS) $(PSBLAS_LIB) $(BLACS) - -srctst: srctst.o - $(MPF90) -o srctst srctst.o $(PSBLAS_LIB) $(BLACS) -.f90.o: - $(MPF90) $(F90COPT) $(INCDIRS) -c $< - -clean: - /bin/rm -f $(DFOBJS) $(ZFOBJS) $(IOOBJS) $(MMHBOBJS) \ - *$(.mod) $(EXEDIR)/df_sample $(EXEDIR)/zf_sample dhb2mm zhb2mm dmm2hb zmm2hb - -lib: - (cd ../../; make library) -verycleanlib: - (cd ../../; make veryclean) - diff --git a/test/Fileread/hbio.f90 b/test/Fileread/hbio.f90 deleted file mode 100644 index 974a2aa9..00000000 --- a/test/Fileread/hbio.f90 +++ /dev/null @@ -1,638 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -module hbio - use psb_sparse_mod - public hb_read, hb_write - interface hb_read - module procedure dhb_read, zhb_read - end interface - interface hb_write - module procedure dhb_write,zhb_write - end interface - -contains - - subroutine dhb_read(a, iret, iunit, filename,b,mtitle) - use psb_sparse_mod - implicit none - type(psb_dspmat_type), intent(out) :: a - integer, intent(out) :: iret - integer, optional, intent(in) :: iunit - character(len=*), optional, intent(in) :: filename - real(kind(1.d0)), optional, pointer :: b(:) - character(len=72), optional, intent(out) :: mtitle - - character :: rhstype,type*3,key*8 - character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20 - integer :: indcrd, ptrcrd, totcrd,& - & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix - integer :: ircode, i,iel,ptr,nzr,infile, j, info - character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)' - character(len=*), parameter :: fmt11='(a1,13x,2i14)' - character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)' - logical, parameter :: debug=.false. - - 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=fmt10) mtitle,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,& - & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt - if (rhscrd.gt.0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix - - call psb_sp_all(a,nnzero,nrow+1,nnzero,ircode) - if (ircode /= 0 ) then - write(0,*) 'Memory allocation failed' - goto 993 - end if - - a%m = nrow - a%k = ncol - a%fida = 'CSR' - a%descra='G' - - - if (tolower(type(1:1)) == 'r') then - if (tolower(type(2:2)) == 'u') then - - - read (infile,fmt=ptrfmt) (a%ia2(i),i=1,nrow+1) - read (infile,fmt=indfmt) (a%ia1(i),i=1,nnzero) - if (valcrd.gt.0) read (infile,fmt=valfmt) (a%aspk(i),i=1,nnzero) - - if (present(b) .and. (rhscrd.gt.0)) then - if (associated(b)) then - if (size(b) < nrow) deallocate(b) - endif - if (.not.associated(b)) then - allocate(b(nrow),stat=info) - endif - read (infile,fmt=rhsfmt) (b(i),i=1,nrow) - endif - - else if (tolower(type(2:2)) == 's') then - - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read - - read (infile,fmt=ptrfmt) (a%ia2(i),i=1,nrow+1) - read (infile,fmt=indfmt) (a%ia1(i),i=1,nnzero) - if (valcrd.gt.0) read (infile,fmt=valfmt) (a%aspk(i),i=1,nnzero) - - if (present(b) .and. (rhscrd.gt.0)) then - if (associated(b)) then - if (size(b) < nrow) deallocate(b) - endif - if (.not.associated(b)) then - allocate(b(nrow),stat=info) - endif - read (infile,fmt=rhsfmt) (b(i),i=1,nrow) - endif - - call psb_ipcsr2coo(a,ircode) - if (ircode /= 0) goto 993 - - call psb_sp_reall(a,2*nnzero,ircode) - ! A is now in COO format - nzr = nnzero - do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then - nzr = nzr + 1 - a%aspk(nzr) = a%aspk(i) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) - end if - end do - a%infoa(psb_nnz_) = nzr - call psb_ipcoo2csr(a,ircode) - if (ircode /= 0) goto 993 - - else - write(0,*) 'read_matrix: matrix type not yet supported' - iret=904 - end if - 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,*) 'DHB_READ: Unexpected end of file ' - return -993 iret=993 - write(0,*) 'DHB_READ: Memory allocation failure' - return - end subroutine dhb_read - - - - - subroutine dhb_write(a,iret,eiout,filename,key,rhs,mtitle) - use psb_sparse_mod - implicit none - type(psb_dspmat_type), intent(in) :: a - integer, intent(out) :: iret - character(len=*), optional, intent(in) :: mtitle - integer, optional, intent(in) :: eiout - character(len=*), optional, intent(in) :: filename - character(len=*), optional, intent(in) :: key - real(kind(1.d0)), optional :: rhs(:) - integer :: iout - - character(len=*), parameter:: ptrfmt='(10I8)',indfmt='(10I8)' - integer, parameter :: jptr=10,jind=10 - character(len=*), parameter:: valfmt='(4E20.12)',rhsfmt='(4E20.12)' - integer, parameter :: jval=4,jrhs=4 - character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)' - character(len=*), parameter :: fmt11='(a1,13x,2i14)' - character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)' - - character(len=72) :: mtitle_ - character(len=8) :: key_ - - character :: rhstype,type*3 - - integer :: i,indcrd,nrhsvl,ptrcrd,rhscrd,totcrd,valcrd,& - & nrow,ncol,nnzero, neltvl, nrhs, nrhsix - - iret = 0 - - if (present(filename)) then - if (filename=='-') then - iout=6 - else - if (present(eiout)) then - iout = eiout - else - iout=99 - endif - open(iout,file=filename, err=901, action='WRITE') - endif - else - if (present(eiout)) then - iout = eiout - else - iout=6 - endif - endif - - if (present(mtitle)) then - mtitle_ = mtitle - else - mtitle_ = 'Temporary PSBLAS title ' - endif - if (present(key)) then - key_ = key - else - key_ = 'PSBMAT00' - endif - - if (toupper(a%fida) == 'CSR') then - - nrow = a%m - ncol = a%k - nnzero = a%ia2(nrow+1)-1 - - neltvl = 0 - - ptrcrd = (nrow+1)/jptr - if (mod(nrow+1,jptr).gt.0) ptrcrd = ptrcrd + 1 - indcrd = nnzero/jind - if (mod(nnzero,jind).gt.0) indcrd = indcrd + 1 - valcrd = nnzero/jval - if (mod(nnzero,jval).gt.0) valcrd = valcrd + 1 - if (present(rhs)) then - if (size(rhs)1 then we have an overlap in the data - ! distribution. - ! - ! integer :: ictxt - ! 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)), optional :: b_glob(:) - ! on entry: this contains right hand side. - ! on exit : - ! - ! real(kind(1.d0)), allocatable, 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 psb_sparse_mod - implicit none - - ! parameters - type(psb_dspmat_type) :: a_glob - real(kind(1.d0)) :: b_glob(:) - integer :: ictxt - type(psb_dspmat_type) :: a - real(kind(1.d0)), allocatable :: b(:) - 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 - integer, intent(in) :: global_indx, n, np - integer, intent(out) :: nv - integer, intent(out) :: pv(*) - end subroutine parts - end interface - - ! local variables - integer :: np, iam - integer :: ircode, length_row, i_count, j_count,& - & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) - integer, allocatable :: iwork(:) - character :: afmt*5, atyp*5 - integer, allocatable :: irow(:),icol(:) - real(kind(1.d0)), allocatable :: val(:) - integer, parameter :: nb=30 - real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5 - logical, parameter :: newt=.true. - 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 psb_info(ictxt, iam, np) - - if (iam == root) then - ! extract information from a_glob - if (a_glob%fida.ne. 'CSR') then - info=135 - ch_err='CSR' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - nrow = a_glob%m - ncol = a_glob%k - if (nrow /= ncol) then - write(0,*) 'a rectangular matrix ? ',nrow,ncol - info=-1 - call psb_errpush(info,name) - goto 9999 - endif - nnzero = size(a_glob%aspk) - nrhs = 1 - ! broadcast informations to other processors - endif - call psb_bcast(ictxt, nrow,root) - call psb_bcast(ictxt, ncol,root) - call psb_bcast(ictxt, nnzero,root) - call psb_bcast(ictxt, nrhs,root) - - liwork = max(np, nrow + ncol) - allocate(iwork(liwork), stat = info) - if (info /= 0) then - info=2025 - int_err(1)=liwork - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif - if (newt) then - call psb_cdall(nrow,nrow,parts,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_cdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_cdall(nrow,nrow,parts,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_pscdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - call psb_spall(a,desc_a,info,nnz=nnzero/np) - if(info/=0) then - info=4010 - ch_err='psb_psspall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geall(b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_psdsall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - isize = max(3*nb,ncol) - - - allocate(val(nb*ncol),irow(nb*ncol),icol(nb*ncol),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - i_count = 1 - - do while (i_count.le.nrow) - - call parts(i_count,nrow,np,iwork, length_row) - - if (length_row.eq.1) then - j_count = i_count - iproc = iwork(1) - do - j_count = j_count + 1 - if (j_count-i_count >= nb) exit - if (j_count > nrow) exit - call parts(j_count,nrow,np,iwork, length_row) - if (length_row /= 1 ) exit - if (iwork(1) /= iproc ) exit - end do - - ! now we should insert rows i_count..j_count-1 - nnr = j_count - i_count - - if (iam == root) then - - ll = 0 - do i= i_count, j_count-1 - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (iproc == iam) then - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geins(nnr,(/(i,i=i_count,j_count-1)/),& - & b_glob(i_count:j_count-1),b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_ins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,nnr,iproc) - call psb_snd(ictxt,ll,iproc) - call psb_snd(ictxt,irow(1:ll),iproc) - call psb_snd(ictxt,icol(1:ll),iproc) - call psb_snd(ictxt,val(1:ll),iproc) - call psb_snd(ictxt,b_glob(i_count:j_count-1),iproc) - call psb_rcv(ictxt,ll,iproc) - endif - else if (iam /= root) then - - if (iproc == iam) then - call psb_rcv(ictxt,nnr,root) - call psb_rcv(ictxt,ll,root) - if (ll > size(irow)) then - write(0,*) iam,'need to reallocate ',ll - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - endif - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count:i_count+nnr-1),root) - call psb_snd(ictxt,ll,root) - call psb_spins(ll,irow,icol,val,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 psb_geins(nnr,(/(i,i=i_count,i_count+nnr-1)/),& - & b_glob(i_count:i_count+nnr-1),b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - - i_count = j_count - - else - write(0,*) iam,'unexpected turn' - ! here processors are counted 1..np - do j_count = 1, length_row - k_count = iwork(j_count) - if (iam == root) then - - ll = 0 - do i= i_count, i_count - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (k_count == iam) then - - call psb_spins(ll,irow,icol,val,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 psb_geins(1,(/i_count/),b_glob(i_count:i_count),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,ll,k_count) - call psb_snd(ictxt,irow(1:ll),k_count) - call psb_snd(ictxt,icol(1:ll),k_count) - call psb_snd(ictxt,val(1:ll),k_count) - call psb_snd(ictxt,b_glob(i_count),k_count) - call psb_rcv(ictxt,ll,k_count) - endif - else if (iam /= root) then - if (k_count == iam) then - - call psb_rcv(ictxt,ll,root) - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count),root) - call psb_snd(ictxt,ll,root) - call psb_spins(ll,irow,icol,val,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 psb_geins(1,(/i_count/),b_glob(i_count:i_count),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - end do - i_count = i_count + 1 - endif - end do - - if (present(fmt)) then - afmt=fmt - else - afmt = 'CSR' - endif - if (newt) then - - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdasb(desc_a,info) - t1 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_cdasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ictxt) - t2 = psb_wtime() - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - t3 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_spasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - if (iam == root) then - write(*,*) 'descriptor assembly: ',t1-t0 - write(*,*) 'sparse matrix assembly: ',t3-t2 - end if - - - else - call psb_spasb(a,desc_a,info,afmt=afmt,dupl=psb_dupl_err_) - if(info/=0)then - info=4010 - ch_err='psspasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - - - call psb_geasb(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 - deallocate(val,irow,icol,stat=info) - if(info/=0)then - info=4010 - ch_err='deallocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(iwork) - if (iam == root) write (*, fmt = *) 'end matdist' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error(ictxt) - return - end if - return - - end subroutine dmatdistf - - - subroutine dmatdistv (a_glob, a, v, ictxt, 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 - ! 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 :: ictxt - ! 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)), optional :: b_glob(:) - ! on entry: this contains right hand side. - ! on exit : - ! - ! real(kind(1.d0)), allocatable, 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 psb_sparse_mod - implicit none ! parameters - type(psb_dspmat_type) :: a_glob - real(kind(1.d0)) :: b_glob(:) - integer :: ictxt, v(:) - type(psb_dspmat_type) :: a - real(kind(1.d0)), allocatable :: b(:) - type (psb_desc_type) :: desc_a - integer, intent(out) :: info - integer, optional :: inroot - character(len=5), optional :: fmt - - integer :: np, iam - integer :: ircode, length_row, i_count, j_count,& - & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) - integer, allocatable :: iwork(:) - character :: afmt*5, atyp*5 - integer, allocatable :: irow(:),icol(:) - real(kind(1.d0)), allocatable :: val(:) - integer, parameter :: nb=30 - logical, parameter :: newt=.true. - real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5 - character(len=20) :: name, ch_err - - info = 0 - err = 0 - name = 'mat_distv' - call psb_erractionsave(err_act) - - ! executable statements - if (present(inroot)) then - root = inroot - else - root = 0 - end if - - call psb_info(ictxt, iam, np) - if (iam == root) then - ! extract information from a_glob - if (a_glob%fida.ne. 'CSR') then - info=135 - ch_err='CSR' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - nrow = a_glob%m - ncol = a_glob%k - if (nrow /= ncol) then - write(0,*) 'a rectangular matrix ? ',nrow,ncol - info=-1 - call psb_errpush(info,name) - goto 9999 - endif - - nnzero = size(a_glob%aspk) - nrhs = 1 - end if - ! broadcast informations to other processors - call psb_bcast(ictxt,nrow, root) - call psb_bcast(ictxt,ncol, root) - call psb_bcast(ictxt,nnzero, root) - call psb_bcast(ictxt,nrhs, root) - liwork = max(np, nrow + ncol) - allocate(iwork(liwork), 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) - goto 9999 - endif - - call psb_cdall(nrow,v,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_cdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_spall(a,desc_a,info,nnz=((nnzero+np-1)/np)) - if(info/=0) then - info=4010 - ch_err='psb_psspall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geall(b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_psdsall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - isize = max(3*nb,ncol) - - - allocate(val(nb*ncol),irow(nb*ncol),icol(nb*ncol),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - i_count = 1 - - do while (i_count <= nrow) - - j_count = i_count - iproc = v(i_count) - - do - j_count = j_count + 1 - if (j_count-i_count >= nb) exit - if (j_count > nrow) exit - if (v(j_count) /= iproc ) exit - end do - - ! now we should insert rows i_count..j_count-1 - nnr = j_count - i_count - - if (iam == root) then - - ll = 0 - do i= i_count, j_count-1 - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (iproc == iam) then - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_geins(nnr,(/(i,i=i_count,j_count-1)/),b_glob(i_count:j_count-1),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='dsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,nnr,iproc) - call psb_snd(ictxt,ll,iproc) - call psb_snd(ictxt,irow(1:ll),iproc) - call psb_snd(ictxt,icol(1:ll),iproc) - call psb_snd(ictxt,val(1:ll),iproc) - call psb_snd(ictxt,b_glob(i_count:j_count-1),iproc) - call psb_rcv(ictxt,ll,iproc) - endif - else if (iam /= root) then - - if (iproc == iam) then - call psb_rcv(ictxt,nnr,root) - call psb_rcv(ictxt,ll,root) - if (ll > size(val)) then - write(0,*) iam,'need to reallocate ',ll - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count:i_count+nnr-1),root) - call psb_snd(ictxt,ll,root) - - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geins(nnr,(/(i,i=i_count,i_count+nnr-1)/),& - & b_glob(i_count:i_count+nnr-1),b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - i_count = j_count - - 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 psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdasb(desc_a,info) - t1 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_cdasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ictxt) - t2 = psb_wtime() - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - t3 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_spasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_geasb(b,desc_a,info) - - if (iam == root) then - write(*,'("Descriptor assembly : ",es10.4)')t1-t0 - write(*,'("Sparse matrix assembly: ",es10.4)')t3-t2 - end if - - if(info/=0)then - info=4010 - ch_err='psdsasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(iwork) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error(ictxt) - return - end if - return - - end subroutine dmatdistv - - - subroutine zmatdistf (a_glob, a, parts, ictxt, 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 - ! 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 :: ictxt - ! 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)), optional :: b_glob(:) - ! on entry: this contains right hand side. - ! on exit : - ! - ! real(kind(1.d0)), allocatable, 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 psb_sparse_mod - implicit none - - ! parameters - type(psb_zspmat_type) :: a_glob - complex(kind(1.d0)) :: b_glob(:) - integer :: ictxt - type(psb_zspmat_type) :: a - complex(kind(1.d0)), allocatable :: b(:) - 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 - integer, intent(in) :: global_indx, n, np - integer, intent(out) :: nv - integer, intent(out) :: pv(*) - end subroutine parts - end interface - - ! local variables - integer :: np, iam - integer :: ircode, length_row, i_count, j_count,& - & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) - integer, allocatable :: iwork(:) - character :: afmt*5, atyp*5 - integer, allocatable :: irow(:),icol(:) - complex(kind(1.d0)), allocatable :: val(:) - integer, parameter :: nb=30 - real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5 - logical, parameter :: newt=.true. - 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 psb_info(ictxt, iam, np) - if (iam == root) then - ! extract information from a_glob - if (a_glob%fida.ne. 'CSR') then - info=135 - ch_err='CSR' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - nrow = a_glob%m - ncol = a_glob%k - if (nrow /= ncol) then - write(0,*) 'a rectangular matrix ? ',nrow,ncol - info=-1 - call psb_errpush(info,name) - goto 9999 - endif - nnzero = size(a_glob%aspk) - nrhs = 1 - endif - ! broadcast informations to other processors - call psb_bcast(ictxt,nrow, root) - call psb_bcast(ictxt,ncol, root) - call psb_bcast(ictxt,nnzero, root) - call psb_bcast(ictxt,nrhs, root) - liwork = max(np, nrow + ncol) - allocate(iwork(liwork), stat = info) - if (info /= 0) then - info=2025 - int_err(1)=liwork - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif - if (newt) then - call psb_cdall(nrow,nrow,parts,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_cdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_cdall(nrow,nrow,parts,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_pscdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - call psb_spall(a,desc_a,info,nnz=nnzero/np) - if(info/=0) then - info=4010 - ch_err='psb_psspall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geall(b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_psdsall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - isize = max(3*nb,ncol) - - - allocate(val(nb*ncol),irow(nb*ncol),icol(nb*ncol),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - i_count = 1 - - do while (i_count.le.nrow) - - call parts(i_count,nrow,np,iwork, length_row) - - if (length_row.eq.1) then - j_count = i_count - iproc = iwork(1) - do - j_count = j_count + 1 - if (j_count-i_count >= nb) exit - if (j_count > nrow) exit - call parts(j_count,nrow,np,iwork, length_row) - if (length_row /= 1 ) exit - if (iwork(1) /= iproc ) exit - end do - - ! now we should insert rows i_count..j_count-1 - nnr = j_count - i_count - - if (iam == root) then - - ll = 0 - do i= i_count, j_count-1 - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (iproc == iam) then - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geins(nnr,(/(i,i=i_count,j_count-1)/),b_glob(i_count:j_count-1),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_ins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,nnr,iproc) - call psb_snd(ictxt,ll,iproc) - call psb_snd(ictxt,irow(1:ll),iproc) - call psb_snd(ictxt,icol(1:ll),iproc) - call psb_snd(ictxt,val(1:ll),iproc) - call psb_snd(ictxt,b_glob(i_count:j_count-1),iproc) - call psb_rcv(ictxt,ll,iproc) - endif - else if (iam /= root) then - - if (iproc == iam) then - call psb_rcv(ictxt,nnr,root) - call psb_rcv(ictxt,ll,root) - if (ll > size(irow)) then - write(0,*) iam,'need to reallocate ',ll - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - endif - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count:i_count+nnr-1),root) - call psb_snd(ictxt,ll,root) - call psb_spins(ll,irow,icol,val,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 psb_geins(nnr,(/(i,i=i_count,i_count+nnr-1)/),& - & b_glob(i_count:i_count+nnr-1),b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - - i_count = j_count - - else - write(0,*) iam,'unexpected turn' - ! here processors are counted 1..np - do j_count = 1, length_row - k_count = iwork(j_count) - if (iam == root) then - - ll = 0 - do i= i_count, i_count - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (k_count == iam) then - - call psb_spins(ll,irow,icol,val,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 psb_geins(1,(/i_count/),b_glob(i_count:i_count),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,ll,k_count) - call psb_snd(ictxt,irow(1:ll),k_count) - call psb_snd(ictxt,icol(1:ll),k_count) - call psb_snd(ictxt,val(1:ll),k_count) - call psb_snd(ictxt,b_glob(i_count),k_count) - call psb_rcv(ictxt,ll,k_count) - endif - else if (iam /= root) then - if (k_count == iam) then - call psb_rcv(ictxt,ll,root) - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count),root) - call psb_snd(ictxt,ll,root) - call psb_spins(ll,irow,icol,val,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 psb_geins(1,(/i_count/),b_glob(i_count:i_count),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - end do - i_count = i_count + 1 - endif - end do - - if (present(fmt)) then - afmt=fmt - else - afmt = 'CSR' - endif - if (newt) then - - call psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdasb(desc_a,info) - t1 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_cdasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ictxt) - t2 = psb_wtime() - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - t3 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_spasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - if (iam == root) then - write(*,*) 'descriptor assembly: ',t1-t0 - write(*,*) 'sparse matrix assembly: ',t3-t2 - end if - - - else - call psb_spasb(a,desc_a,info,afmt=afmt,dupl=psb_dupl_err_) - if(info/=0)then - info=4010 - ch_err='psspasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - - - call psb_geasb(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 - deallocate(val,irow,icol,stat=info) - if(info/=0)then - info=4010 - ch_err='deallocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(iwork) - if (iam == root) write (*, fmt = *) 'end matdist' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error(ictxt) - return - end if - return - - end subroutine zmatdistf - - - subroutine zmatdistv (a_glob, a, v, ictxt, 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 - ! 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 :: ictxt - ! 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)), optional :: b_glob(:) - ! on entry: this contains right hand side. - ! on exit : - ! - ! real(kind(1.d0)), allocatable, 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 psb_sparse_mod - implicit none ! parameters - type(psb_zspmat_type) :: a_glob - complex(kind(1.d0)) :: b_glob(:) - integer :: ictxt, v(:) - type(psb_zspmat_type) :: a - complex(kind(1.d0)), allocatable :: b(:) - type(psb_desc_type) :: desc_a - integer, intent(out) :: info - integer, optional :: inroot - character(len=5), optional :: fmt - - integer :: np, iam - integer :: ircode, length_row, i_count, j_count,& - & k_count, blockdim, root, liwork, nrow, ncol, nnzero, nrhs,& - & i,j,k, ll, nz, isize, iproc, nnr, err, err_act, int_err(5) - integer, allocatable :: iwork(:) - character :: afmt*5, atyp*5 - integer, allocatable :: irow(:),icol(:) - complex(kind(1.d0)), allocatable :: val(:) - integer, parameter :: nb=30 - logical, parameter :: newt=.true. - real(kind(1.d0)) :: t0, t1, t2, t3, t4, t5 - character(len=20) :: name, ch_err - - info = 0 - err = 0 - name = 'mat_distv' - call psb_erractionsave(err_act) - - ! executable statements - if (present(inroot)) then - root = inroot - else - root = 0 - end if - - call psb_info(ictxt, iam, np) - if (iam == root) then - ! extract information from a_glob - if (toupper(a_glob%fida) /= 'CSR') then - info=135 - ch_err='CSR' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - nrow = a_glob%m - ncol = a_glob%k - if (nrow /= ncol) then - write(0,*) 'a rectangular matrix ? ',nrow,ncol - info=-1 - call psb_errpush(info,name) - goto 9999 - endif - - nnzero = size(a_glob%aspk) - nrhs = 1 - end if - ! broadcast informations to other processors - call psb_bcast(ictxt,nrow, root) - call psb_bcast(ictxt,ncol, root) - call psb_bcast(ictxt,nnzero, root) - call psb_bcast(ictxt,nrhs, root) - liwork = max(np, nrow + ncol) - allocate(iwork(liwork), 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) - goto 9999 - endif - - call psb_cdall(nrow,v,ictxt,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_cdall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_spall(a,desc_a,info,nnz=((nnzero+np-1)/np)) - if(info/=0) then - info=4010 - ch_err='psb_psspall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geall(b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_psdsall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - isize = max(3*nb,ncol) - - - allocate(val(nb*ncol),irow(nb*ncol),icol(nb*ncol),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - i_count = 1 - - do while (i_count <= nrow) - - j_count = i_count - iproc = v(i_count) - - do - j_count = j_count + 1 - if (j_count-i_count >= nb) exit - if (j_count > nrow) exit - if (v(j_count) /= iproc ) exit - end do - - ! now we should insert rows i_count..j_count-1 - nnr = j_count - i_count - - if (iam == root) then - ll = 0 - do i= i_count, j_count-1 - call psb_sp_getrow(i,a_glob,nz,& - & irow(ll+1:),icol(ll+1:),val(ll+1:), info) - if (info /= 0) then - if (nz >min(size(irow(ll+1:)),size(icol(ll+1:)),size(val(ll+1:)))) then - write(0,*) 'Allocation failure? This should not happen!' - end if - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ll = ll + nz - end do - - if (iproc == iam) then - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psb_spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_geins(nnr,(/(i,i=i_count,j_count-1)/),b_glob(i_count:j_count-1),& - & b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='dsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - else - call psb_snd(ictxt,nnr,iproc) - call psb_snd(ictxt,ll,iproc) - call psb_snd(ictxt,irow(1:ll),iproc) - call psb_snd(ictxt,icol(1:ll),iproc) - call psb_snd(ictxt,val(1:ll),iproc) - call psb_snd(ictxt,b_glob(i_count:j_count-1),iproc) - call psb_rcv(ictxt,ll,iproc) - endif - else if (iam /= root) then - - if (iproc == iam) then - call psb_rcv(ictxt,nnr,root) - call psb_rcv(ictxt,ll,root) - if (ll > size(val)) then - write(0,*) iam,'need to reallocate ',ll - deallocate(val,irow,icol) - allocate(val(ll),irow(ll),icol(ll),stat=info) - if(info/=0) then - info=4010 - ch_err='Allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - call psb_rcv(ictxt,irow(1:ll),root) - call psb_rcv(ictxt,icol(1:ll),root) - call psb_rcv(ictxt,val(1:ll),root) - call psb_rcv(ictxt,b_glob(i_count:i_count+nnr-1),root) - call psb_snd(ictxt,ll,root) - - call psb_spins(ll,irow,icol,val,a,desc_a,info) - if(info/=0) then - info=4010 - ch_err='spins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - call psb_geins(nnr,(/(i,i=i_count,i_count+nnr-1)/),& - & b_glob(i_count:i_count+nnr-1),b,desc_a,info) - if(info/=0) then - info=4010 - ch_err='psdsins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - endif - i_count = j_count - - 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 psb_barrier(ictxt) - t0 = psb_wtime() - call psb_cdasb(desc_a,info) - t1 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_cdasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_barrier(ictxt) - t2 = psb_wtime() - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - t3 = psb_wtime() - if(info/=0)then - info=4010 - ch_err='psb_spasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_geasb(b,desc_a,info) - - if (iam == root) then - write(*,'("Descriptor assembly : ",es10.4)')t1-t0 - write(*,'("Sparse matrix assembly: ",es10.4)')t3-t2 - end if - - if(info/=0)then - info=4010 - ch_err='psdsasb' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(iwork) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error(ictxt) - return - end if - return - - end subroutine zmatdistv - -end module mat_dist diff --git a/test/Fileread/mmio.f90 b/test/Fileread/mmio.f90 deleted file mode 100644 index da034346..00000000 --- a/test/Fileread/mmio.f90 +++ /dev/null @@ -1,379 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -module mmio - use psb_sparse_mod - public mm_mat_read, mm_mat_write - interface mm_mat_read - module procedure dmm_mat_read, zmm_mat_read - end interface - interface mm_mat_write - module procedure dmm_mat_write,zmm_mat_write - end interface - -contains - - subroutine dmm_mat_read(a, iret, iunit, filename) - use psb_sparse_mod - implicit none - type(psb_dspmat_type), 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 - character(1024) :: line - integer :: nrow, ncol, nnzero, neltvl, nrhs, nrhsix - integer :: ircode, i,iel,nzr,infile, j - logical, parameter :: debug=.false. - - 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 - - if ( (tolower(object) /= 'matrix').or.(tolower(fmt)/='coordinate')) then - write(0,*) 'READ_MATRIX: input file type not yet supported' - iret=909 - return - end if - if (debug) write(*,*) mmheader,':', object, ':',fmt,':', type,':', sym - - do - read(infile,fmt='(a)') line - if (line(1:1) /= '%') exit - end do - if (debug) write(*,*) 'Line on input : "',line,'"' - read(line,fmt=*) nrow,ncol,nnzero - if (debug) write(*,*) 'Out: ',nrow,ncol,nnzero - - if ((tolower(type) == 'real').and.(tolower(sym) == 'general')) then - call psb_sp_all(nrow,ncol,a,nnzero,ircode) - a%fida = 'COO' - a%descra = 'G' - 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 - a%infoa(psb_nnz_) = nnzero - call psb_ipcoo2csr(a,ircode) - - else if ((tolower(type) == 'real').and.(tolower(sym) == 'symmetric')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read - call psb_sp_all(nrow,ncol,a,2*nnzero,ircode) - a%fida = 'COO' - a%descra = 'G' - if (ircode /= 0) goto 993 - do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),a%aspk(i) - end do - - nzr = nnzero - do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then - nzr = nzr + 1 - a%aspk(nzr) = a%aspk(i) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) - end if - end do - a%infoa(psb_nnz_) = nzr - call psb_ipcoo2csr(a,ircode) - - 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 dmm_mat_read - - - - - - subroutine dmm_mat_write(a,mtitle,iret,eiout,filename) - use psb_sparse_mod - implicit none - type(psb_dspmat_type), intent(in) :: a - integer, intent(out) :: iret - character(len=*), intent(in) :: mtitle - integer, optional, intent(in) :: eiout - character(len=*), optional, intent(in) :: filename - integer :: iout - - - iret = 0 - - if (present(filename)) then - if (filename=='-') then - iout=6 - else - if (present(eiout)) then - iout = eiout - else - iout=99 - endif - open(iout,file=filename, err=901, action='WRITE') - endif - else - if (present(eiout)) then - iout = eiout - else - iout=6 - endif - endif - - call psb_csprt(iout,a,head=mtitle) - - if (iout /= 6) close(iout) - - - return - -901 continue - iret=901 - write(0,*) 'Error while opening ',filename - return - end subroutine dmm_mat_write - - - - subroutine zmm_mat_read(a, iret, iunit, filename) - use psb_sparse_mod - implicit none - type(psb_zspmat_type), 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 - character(1024) :: line - integer :: nrow, ncol, nnzero, neltvl, nrhs, nrhsix - integer :: ircode, i,iel,nzr,infile,j - real(kind(1.d0)) :: are, aim - logical, parameter :: debug=.false. - - - 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 - - if ( (tolower(object) /= 'matrix').or.(tolower(fmt)/='coordinate')) then - write(0,*) 'READ_MATRIX: input file type not yet supported' - iret=909 - return - end if - if (debug) write(*,*) mmheader,':', object, ':',fmt,':', type,':', sym - - do - read(infile,fmt='(a)') line - if (line(1:1) /= '%') exit - end do - if (debug) write(*,*) 'Line on input : "',line,'"' - read(line,fmt=*) nrow,ncol,nnzero - if (debug) write(*,*) 'Out: ',nrow,ncol,nnzero - - if ((tolower(type) == 'complex').and.(tolower(sym) == 'general')) then - call psb_sp_all(nrow,ncol,a,nnzero,ircode) - if (ircode /= 0) goto 993 - a%fida = 'COO' - a%descra = 'G' - do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim - a%aspk(i) = cmplx(are,aim) - end do - a%infoa(psb_nnz_) = nnzero - - call psb_ipcoo2csr(a,ircode) - - else if ((tolower(type) == 'complex').and.(tolower(sym) == 'symmetric')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read - call psb_sp_all(nrow,ncol,a,2*nnzero,ircode) - if (ircode /= 0) goto 993 - a%fida = 'COO' - a%descra = 'G' - do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim - a%aspk(i) = cmplx(are,aim) - end do - - nzr = nnzero - do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then - nzr = nzr + 1 - a%aspk(nzr) = a%aspk(i) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) - end if - end do - a%infoa(psb_nnz_) = nzr - call psb_ipcoo2csr(a,ircode) - - else if ((tolower(type) == 'complex').and.(tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read - call psb_sp_all(nrow,ncol,a,2*nnzero,ircode) - if (ircode /= 0) goto 993 - a%fida = 'COO' - a%descra = 'G' - do i=1,nnzero - read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim - a%aspk(i) = cmplx(are,aim) - end do - - nzr = nnzero - do i=1,nnzero - if (a%ia1(i) /= a%ia2(i)) then - nzr = nzr + 1 - a%aspk(nzr) = conjg(a%aspk(i)) - a%ia1(nzr) = a%ia2(i) - a%ia2(nzr) = a%ia1(i) - end if - end do - a%infoa(psb_nnz_) = nzr - call psb_ipcoo2csr(a,ircode) - - 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 zmm_mat_write(a,mtitle,iret,eiout,filename) - use psb_sparse_mod - implicit none - type(psb_zspmat_type), intent(in) :: a - integer, intent(out) :: iret - character(len=*), intent(in) :: mtitle - integer, optional, intent(in) :: eiout - character(len=*), optional, intent(in) :: filename - integer :: iout - - - iret = 0 - - if (present(filename)) then - if (filename=='-') then - iout=6 - else - if (present(eiout)) then - iout = eiout - else - iout=99 - endif - open(iout,file=filename, err=901, action='WRITE') - endif - else - if (present(eiout)) then - iout = eiout - else - iout=6 - endif - endif - - call psb_csprt(iout,a,head=mtitle) - - if (iout /= 6) close(iout) - - - return - -901 continue - iret=901 - write(0,*) 'Error while opening ',filename - return - end subroutine zmm_mat_write - - -end module mmio diff --git a/test/Fileread/part_blk2.f b/test/Fileread/part_blk2.f deleted file mode 100644 index 22ef34f5..00000000 --- a/test/Fileread/part_blk2.f +++ /dev/null @@ -1,86 +0,0 @@ -C -C Parallel Sparse BLAS v2.0 -C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -C Alfredo Buttari University of Rome Tor Vergata -C -C Redistribution and use in source and binary forms, with or without -C modification, are permitted provided that the following conditions -C are met: -C 1. Redistributions of source code must retain the above copyright -C notice, this list of conditions and the following disclaimer. -C 2. Redistributions in binary form must reproduce the above copyright -C notice, this list of conditions, and the following disclaimer in the -C documentation and/or other materials provided with the distribution. -C 3. The name of the PSBLAS group or the names of its contributors may -C not be used to endorse or promote products derived from this -C software without specific written permission. -C -C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -C POSSIBILITY OF SUCH DAMAGE. -C -C -C -C User defined function corresponding to an HPF BLOCK partition -C - SUBROUTINE PART_BLK2(IDX,N,NP,PV,NV) - - IMPLICIT NONE - - INTEGER IDX, N, NP - INTEGER NV - INTEGER PV(*) - DOUBLE PRECISION DDIFF - INTEGER IB1, IB2, IP, NB, NB1, NNB1 - - NV = 1 - NB = N/NP - NB1 = NB+1 - NNB1 = MOD(N,NP) - IF (IDX .LE. (NNB1*NB1)) THEN - PV(1) = (IDX - 1) / NB1 - ELSE - IF (NB > 0) THEN - IP = ( (IDX-NNB1*NB1) - 1)/NB - PV(1) = NNB1 + IP - ELSE - write(0,*) 'Impossible ??? ' - PV(1) = NNB1 - ENDIF - ENDIF - - RETURN - END - - - SUBROUTINE BLD_PARTBLK2(N,NP,IVG) - - INTEGER N, IVG(*),NP - INTEGER IB1, IB2, IP, NB, NB1, NNB1, I - - NB = N/NP - NB1 = NB+1 - NNB1 = MOD(N,NP) - DO I=1,N - IF (I .LE. (NNB1*NB1)) THEN - IVG(I) = (I - 1) / NB1 - ELSE - IF (NB > 0) THEN - IP = ( (I-NNB1*NB1) - 1)/NB - IVG(I) = NNB1 + IP - ELSE - write(0,*) 'Impossible ??? ' - IVG(I) = NNB1 - ENDIF - ENDIF - ENDDO - - END diff --git a/test/Fileread/part_block.f b/test/Fileread/part_block.f deleted file mode 100644 index a4142756..00000000 --- a/test/Fileread/part_block.f +++ /dev/null @@ -1,97 +0,0 @@ -C -C Parallel Sparse BLAS v2.0 -C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -C Alfredo Buttari University of Rome Tor Vergata -C -C Redistribution and use in source and binary forms, with or without -C modification, are permitted provided that the following conditions -C are met: -C 1. Redistributions of source code must retain the above copyright -C notice, this list of conditions and the following disclaimer. -C 2. Redistributions in binary form must reproduce the above copyright -C notice, this list of conditions, and the following disclaimer in the -C documentation and/or other materials provided with the distribution. -C 3. The name of the PSBLAS group or the names of its contributors may -C not be used to endorse or promote products derived from this -C software without specific written permission. -C -C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -C POSSIBILITY OF SUCH DAMAGE. -C -C -C -C User defined function corresponding to an HPF BLOCK partition -C - SUBROUTINE PART_BLOCK(GLOBAL_INDX,N,NP,PV,NV) - - IMPLICIT NONE - - INTEGER GLOBAL_INDX, N, NP - INTEGER NV - INTEGER PV(*) - INTEGER DIM_BLOCK - DOUBLE PRECISION DDIFF - INTEGER IB1, IB2, IPV - - double precision PC - PARAMETER (PC=0.0D0) - - DIM_BLOCK = (N + NP - 1)/NP - NV = 1 - PV(NV) = (GLOBAL_INDX - 1) / DIM_BLOCK - - IPV = PV(1) - IB1 = IPV * DIM_BLOCK + 1 - IB2 = (IPV+1) * DIM_BLOCK - - DDIFF = DBLE(ABS(GLOBAL_INDX-IB1))/DBLE(DIM_BLOCK) - IF (DDIFF .lt. PC/2) THEN -C -C Overlap at the beginning of a block, with the previous proc -C - IF (IPV.gt.0) THEN - NV = NV + 1 - PV(NV) = IPV - 1 - ENDIF - ENDIF - - DDIFF = DBLE(ABS(GLOBAL_INDX-IB2))/DBLE(DIM_BLOCK) - IF (DDIFF .lt. PC/2) THEN -C -C Overlap at the end of a block, with the next proc -C - IF (IPV.lt.(NP-1)) THEN - NV = NV + 1 - PV(NV) = IPV + 1 - ENDIF - ENDIF - - RETURN - END - - - - SUBROUTINE BLD_PARTBLOCK(N,NP,IVG) - - INTEGER N,NP,IVG(*) - - INTEGER DIM_BLOCK,I - - - DIM_BLOCK = (N + NP - 1)/NP - DO I=1,N - IVG(I) = (I - 1) / DIM_BLOCK - ENDDO - - END - - diff --git a/test/Fileread/partgraph.f90 b/test/Fileread/partgraph.f90 deleted file mode 100644 index 6622d8b2..00000000 --- a/test/Fileread/partgraph.f90 +++ /dev/null @@ -1,222 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -! -! Purpose: -! Provide a set of subroutines to define a data distribution based on -! a graph partitioning routine. -! -! Subroutines: -! -! BUILD_GRPPART(A,NPARTS): This subroutine will be called by the root -! process to build define the data distribuition mapping. -! Input parameters: -! TYPE(D_SPMAT) :: A The input matrix. The coefficients are -! ignored; only the structure is used. -! INTEGER :: NPARTS How many parts we are requiring to the -! partition utility -! -! DISTR_GRPPART(RROOT,CROOT,ICTXT): This subroutine will be called by -! all processes to distribute the information computed by the root -! process, to be used subsequently. -! -! -! PART_GRAPH : The subroutine to be passed to PSBLAS sparse library; -! uses information prepared by the previous two subroutines. -! -MODULE PARTGRAPH - public part_graph, build_grppart, distr_grppart,& - & getv_grppart, build_usrpart, free_part - private - integer, allocatable, save :: graph_vect(:) - -CONTAINS - - subroutine part_graph(global_indx,n,np,pv,nv) - - integer, intent(in) :: global_indx, n, np - integer, intent(out) :: nv - integer, intent(out) :: pv(*) - - IF (.not.allocated(graph_vect)) then - write(0,*) 'Fatal error in PART_GRAPH: vector GRAPH_VECT ',& - & 'not initialized' - return - endif - if ((global_indx<1).or.(global_indx > size(graph_vect))) then - write(0,*) 'Fatal error in PART_GRAPH: index GLOBAL_INDX ',& - & 'outside GRAPH_VECT bounds',global_indx,size(graph_vect) - return - endif - nv = 1 - pv(nv) = graph_vect(global_indx) - return - end subroutine part_graph - - - subroutine distr_grppart(root, ictxt) - use psb_sparse_mod - integer :: root, ictxt - integer :: n, me, np - - call psb_info(ictxt,me,np) - - if (.not.((root>=0).and.(root