diff --git a/Make.inc.in b/Make.inc.in index 5ff454bf..7be1d4c5 100755 --- a/Make.inc.in +++ b/Make.inc.in @@ -60,6 +60,7 @@ BASELIBNAME=@BASELIBNAME@ PRECLIBNAME=@PRECLIBNAME@ METHDLIBNAME=@METHDLIBNAME@ UTILLIBNAME=@UTILLIBNAME@ +EIGENLIBNAME=@EIGENLIBNAME@ BASEMODNAME=@BASEMODNAME@ PRECMODNAME=@PRECMODNAME@ METHDMODNAME=@METHDMODNAME@ diff --git a/Makefile b/Makefile index 5479fa6b..699a1642 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ include Make.inc -all: libd based precd kryld utild +all: libd based precd kryld utild eigend @echo "=====================================" @echo "PSBLAS libraries Compilation Successful." @@ -19,6 +19,8 @@ kryld: cd krylov && $(MAKE) lib utild: cd util&& $(MAKE) lib +eigend: + cd eigen&& $(MAKE) lib install: all (./mkdir.sh $(INSTALL_INCLUDEDIR) &&\ @@ -36,6 +38,7 @@ clean: cd prec && $(MAKE) clean cd krylov && $(MAKE) clean cd util && $(MAKE) clean + cd eigen && $(MAKE) clean check: all make check -C test/serial diff --git a/configure b/configure index e3a2dfbb..25436109 100755 --- a/configure +++ b/configure @@ -641,6 +641,7 @@ ac_subst_vars='am__EXEEXT_FALSE am__EXEEXT_TRUE LTLIBOBJS LIBOBJS +EIGENLIBNAME UTILLIBNAME METHDLIBNAME PRECLIBNAME @@ -10933,6 +10934,7 @@ BASELIBNAME=libpsb_base.a PRECLIBNAME=libpsb_prec.a METHDLIBNAME=libpsb_krylov.a UTILLIBNAME=libpsb_util.a +EIGENLIBNAME=libpsb_eigen.a ############################################################################### # Variable substitutions : the Make.inc.in will have these @VARIABLES@ @@ -11010,6 +11012,7 @@ FDEFINES=$(PSBFDEFINES) + ############################################################################### # the following files will be created by Automake diff --git a/configure.ac b/configure.ac index f7771132..5ddcd186 100755 --- a/configure.ac +++ b/configure.ac @@ -740,6 +740,7 @@ BASELIBNAME=libpsb_base.a PRECLIBNAME=libpsb_prec.a METHDLIBNAME=libpsb_krylov.a UTILLIBNAME=libpsb_util.a +EIGENLIBNAME=libpsb_eigen.a ############################################################################### # Variable substitutions : the Make.inc.in will have these @VARIABLES@ @@ -816,6 +817,7 @@ AC_SUBST(BASELIBNAME) AC_SUBST(PRECLIBNAME) AC_SUBST(METHDLIBNAME) AC_SUBST(UTILLIBNAME) +AC_SUBST(EIGENLIBNAME) ############################################################################### # the following files will be created by Automake diff --git a/eigen/Makefile b/eigen/Makefile new file mode 100644 index 00000000..1ef26ca0 --- /dev/null +++ b/eigen/Makefile @@ -0,0 +1,34 @@ +include ../Make.inc + + +HERE=. +LIBDIR=../lib +INCDIR=../include + +MODOBJS=psb_eigen_mod.o +F90OBJS=psb_d_power_vect.o +OBJS=$(F90OBJS) $(MODOBJS) + +LOCAL_MODS=$(MODOBJS:.o=$(.mod)) + +LIBNAME=$(EIGENLIBNAME) + +FINCLUDES=$(FMFLAG). $(FMFLAG)$(INCDIR) + +lib: $(OBJS) + $(AR) $(HERE)/$(LIBNAME) $(OBJS) + $(RANLIB) $(HERE)/$(LIBNAME) + /bin/cp -p $(CPUPDFLAG) $(HERE)/$(LIBNAME) $(LIBDIR) + /bin/cp -p $(CPUPDFLAG) *$(.mod) $(INCDIR) + + +$(F90OBJS): $(MODOBJS) +$(OBJS): $(INCDIR)/$(PRECMODNAME)$(.mod) $(INCDIR)/$(BASEMODNAME)$(.mod) + + +veryclean: clean + /bin/rm -f $(HERE)/$(LIBNAME) + +clean: + /bin/rm -f $(OBJS) *$(.mod) + diff --git a/eigen/psb_d_power_vect.f90 b/eigen/psb_d_power_vect.f90 new file mode 100644 index 00000000..060c57a0 --- /dev/null +++ b/eigen/psb_d_power_vect.f90 @@ -0,0 +1,25 @@ +Subroutine psb_d_power_vect(method,a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,irst,istop,cond) + + use psb_base_mod + use psb_prec_mod + use psb_krylov_mod + use psb_eigen_mod, psb_protect_name =>psb_d_power_vect + + character(len=*) :: method + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(inout) :: prec + type(psb_d_vect_type), Intent(inout) :: b + type(psb_d_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err,cond + + + + +end Subroutine psb_d_power_vect + diff --git a/eigen/psb_eigen_mod.f90 b/eigen/psb_eigen_mod.f90 new file mode 100644 index 00000000..2551b31c --- /dev/null +++ b/eigen/psb_eigen_mod.f90 @@ -0,0 +1,29 @@ +module psb_eigen_mod + use psb_base_mod + use psb_krylov_mod + use psb_prec_mod + + interface + Subroutine psb_d_power_vect(method,a,prec,b,x,eps,desc_a,info,& + & itmax,iter,err,itrace,irst,istop,cond) + + use psb_base_mod, only : psb_ipk_, psb_desc_type, psb_dspmat_type, & + & psb_dpk_, psb_d_vect_type + use psb_prec_mod, only : psb_dprec_type + + character(len=*) :: method + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(psb_dprec_type), intent(inout) :: prec + type(psb_d_vect_type), Intent(inout) :: b + type(psb_d_vect_type), Intent(inout) :: x + Real(psb_dpk_), Intent(in) :: eps + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, irst,istop + integer(psb_ipk_), Optional, Intent(out) :: iter + Real(psb_dpk_), Optional, Intent(out) :: err,cond + + end Subroutine psb_d_power_vect + end interface + +end module psb_eigen_mod diff --git a/test/eigen/Makefile b/test/eigen/Makefile new file mode 100644 index 00000000..314f63e5 --- /dev/null +++ b/test/eigen/Makefile @@ -0,0 +1,84 @@ +# +# Libraries used +# +BASEDIR=../.. +INCDIR=$(BASEDIR)/include/ +include $(INCDIR)/Make.inc.psblas +LIBDIR=$(BASEDIR)/lib/ +PSBLAS_LIB= -L$(LIBDIR) -leigen -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base +LDLIBS=$(PSBLDLIBS) + +FINCLUDES=$(FMFLAG)$(INCDIR) $(FMFLAG). + +DTOBJS=d_file_spmv.o +STOBJS=s_file_spmv.o +DPGOBJS=pdgenspmv.o +PWFOBJS = power_file.o +PWFROBJS = power_file_real.o #eig_mod.o +ARNOBJS =arnoldi_file.o +ARNROBJS =arnoldi_file_real.o +SIROBJS=shift_invert_real.o +EXEDIR=./runs + +all: d_file_spmv s_file_spmv pdgenspmv power_file power_file_real + + +d_file_spmv: $(DTOBJS) + $(F90LINK) $(LOPT) $(DTOBJS) -o d_file_spmv $(PSBLAS_LIB) $(LDLIBS) + /bin/mv d_file_spmv $(EXEDIR) + +pdgenspmv: $(DPGOBJS) + $(F90LINK) $(LOPT) $(DPGOBJS) -o pdgenspmv $(PSBLAS_LIB) $(LDLIBS) + /bin/mv pdgenspmv $(EXEDIR) + +s_file_spmv: $(STOBJS) + $(F90LINK) $(LOPT) $(STOBJS) -o s_file_spmv $(PSBLAS_LIB) $(LDLIBS) + /bin/mv s_file_spmv $(EXEDIR) + +power_file: $(PWFOBJS) + $(F90LINK) $(LOPT) $(PWFOBJS) -o power_file $(PSBLAS_LIB) $(LDLIBS) + /bin/mv power_file $(EXEDIR) + +power_file_real: $(PWFROBJS) + $(F90LINK) $(LOPT) $(PWFROBJS) -o power_file_real $(PSBLAS_LIB) $(LDLIBS) + /bin/mv power_file_real $(EXEDIR) + +arnoldi_file: $(ARNOBJS) + $(F90LINK) $(LOPT) $(ARNOBJS) -o arnoldi_file $(PSBLAS_LIB) $(LDLIBS) + /bin/mv arnoldi_file $(EXEDIR) + +arnoldi_file_real: $(ARNROBJS) . + $(F90LINK) $(LOPT) $(ARNROBJS) -o arnoldi_file_real $(PSBLAS_LIB) $(LDLIBS) + /bin/mv arnoldi_file_real $(EXEDIR) + +shift_invert_real: $(SIROBJS) + $(F90LINK) $(LOPT) $(SIROBJS) -o shift_invert_real $(PSBLAS_LIB) $(LDLIBS) + /bin/mv shift_invert_real $(EXEDIR) + +analyse: analyse.o + $(F90LINK) $(LOPT) analyse.o -o analyse $(PSBLAS_LIB) $(LDLIBS) + /bin/mv analyse /server/tamestoy/TEMP/MaxEigenAdj/ + +an_lapl: an_lapl.o + $(F90LINK) $(LOPT) an_lapl.o -o an_lapl $(PSBLAS_LIB) $(LDLIBS) + /bin/mv an_lapl /server/tamestoy/TEMP/SCALING-BA + +laplacian: laplacian.o + $(F90LINK) $(LOPT) laplacian.o -o laplacian $(PSBLAS_LIB) $(LDLIBS) + /bin/mv laplacian $(EXEDIR) + +test_chseqr: test_chseqr.o + $(F90LINK) $(LOPT) test_chseqr.o -o test_chseqr $(PSBLAS_LIB) $(LDLIBS) + +#eig_mod.o: eig_mod.f90 + #$(F90LINK) $(LOPT) -c eig_mod.f90 $(PSBLAS_LIB) $(LDLIBS) + #gfortran -c eig_mod.f90 + +clean: + /bin/rm -f *.o + #/bin/rm -f $(DBOBJSS) $(DBOBJS) $(PWOBJS) $(PWFOBJS) $(DTOBJS) $(STOBJS) $(ARNOBJS) + +lib: + (cd ../../; make library) +verycleanlib: + (cd ../../; make veryclean) diff --git a/test/eigen/adj_read.f90 b/test/eigen/adj_read.f90 new file mode 100644 index 00000000..853aeee5 --- /dev/null +++ b/test/eigen/adj_read.f90 @@ -0,0 +1,38 @@ +subroutine adj_read (a,filename,desc_a,info) + + type(psb_dspmat_type), intent (out) :: a + character(len=20) :: filename + type(psb_desc_type):: desc_a + integer (psb_ipk_) :: info + implicit none + + integer :: i,j,nnzero,nbCol + integer :: unitFile, iError,line + integer(psb_ipk_), allocatable :: fileContent(:,:),ia(:),ja(:) + real (psb_dpk_), allocatable :: val(:) + integer(psb_ipk_), allocatable :: size_mat(:) + + nbCol = 2 + allocate (size_mat(nbCol)) + unitFile = 1 + open(UNIT=unitFile, FILE=filename, FORM="FORMATTED", STATUS="OLD", +ACTION="READ") + + nnzero = size_mat(2) + allocate (fileContent(nnzeros,nbCol)) + + do line = 1,nnzero + read(unitFile, *) fileContent(line,1:nbCol) + end do saveNodes + close(UNIT=unitFile) + + allocate(ia(nnzero),ja(nnzero),val(nnzero)) + do i=1,nnzero + ia(i)=fileContent(i,1) + ja(i)=filecontent(i,2) + val(i)=1.0 + end do + + call psb_spins(nnzero, ia, ja, val, a, desc_a, info) + +end subroutine adj_read diff --git a/test/eigen/an_lapl.f90 b/test/eigen/an_lapl.f90 new file mode 100644 index 00000000..73a0181e --- /dev/null +++ b/test/eigen/an_lapl.f90 @@ -0,0 +1,82 @@ +program analyse + use psb_base_mod + use psb_util_mod + implicit none + + character(len=30) :: lapl_file,max_file,res_file,analyse_file,mat_file + integer (psb_ipk_) :: i,nb_mat, n + real(psb_dpk_) ::tab(18),eig1,eig2,eig0,eigmin,eigmax,delta,delta2,& + &deltamax,timemin,timemax,t1,t2,time0,t0 + + read(psb_inp_unit,*) lapl_file + read(psb_inp_unit,*) max_file + read(psb_inp_unit,*) res_file + read(psb_inp_unit,*) analyse_file + read(psb_inp_unit,*) nb_mat + + open(15,FILE=lapl_file,action="read") + open(13,FILE=max_file,action="read") + open(14,file=res_file,action="read") + + do i=1,21 + read(15,*) + end do + do i=1,22 + read(13,*) + end do + + open (16, file=analyse_file,position="append",action="write") + n=0 + delta=0 + delta2=0 + deltamax=0 + timemax=0 + timemin=0 + time0=0 + do i=1,nb_mat + read(15,*)tab(1:18) + read(14,*)mat_file,eig1,eig2,t1 + read(14,*)eigmin,eig0,t2 + + if(abs(tab(6)-eigmin)>0.001)then + n=n+1 + write(16,'("smallest of : ",F20.2,F20.2)')tab(2),tab(3) + write(16,'(F20.7,F20.7)')tab(6),eigmin + end if + if(abs(tab(16)-eig2)>0.01)then + n=n+1 + write(16,'("lambda N-1 of : ",F20.2,F20.2)')tab(2),tab(3) + write(16,'(F20.7,F20.7)')tab(16),eig2 + end if + if(abs(tab(18)-eig1)>0.01)then + n=n+1 + write(16,'("biggest of : ", F20.2,F20.2)')tab(2),tab(3) + write(16,'(F20.7,F20.7)')tab(18),eig1 + end if + timemax=timemax+t1 + timemin=timemin+t2 + delta2=delta2+abs(1-eigmin/tab(6)) + delta=delta+abs(1-eig2/tab(16))+abs(1-eig1/tab(18)) + + read(13,*)tab(1:7) + read(14,*)mat_file,eigmax,t0 + if(abs(tab(4)-eigmax)>0.01)then + write(16,'(F20.2,F20.2)')tab(2),tab(3) + write(16,'(F20.6,F20.6)')tab(4),eigmax + endif + time0=time0+t0 + deltamax=deltamax+abs(1-eigmax/tab(4)) + end do + write(16,'("errors on ",i20," eigenvalues")')n + write(16,'("gap average for lapl max",g20.5)')delta/(2*nb_mat) + write(16,'("gap average for lapl min",g20.5)')delta2/(nb_mat) + write(16,'("gap average for max",g20.5)')deltamax/(nb_mat) + write(16,'("time average lapl max",g20.5)')timemax/(nb_mat) + write(16,'("time average lapl min",g20.5)')timemin/(nb_mat) + write(16,'("time average A max",g20.5)')time0/(nb_mat) + close(15) + close(14) + close(13) + close(16) + +end program analyse diff --git a/test/eigen/analyse.f90 b/test/eigen/analyse.f90 new file mode 100644 index 00000000..5dc60edd --- /dev/null +++ b/test/eigen/analyse.f90 @@ -0,0 +1,43 @@ +program analyse + use psb_base_mod + use psb_util_mod + implicit none + + character(len=30) :: given_file,res_file,analyse_file,mat_file + integer (psb_ipk_) :: i,nb_mat, n + real(psb_dpk_) :: tab(7),eig + + read(psb_inp_unit,*) given_file + read(psb_inp_unit,*) res_file + read(psb_inp_unit,*) analyse_file + read(psb_inp_unit,*) nb_mat + + + open(15,FILE=given_file,action="read") + open(14,file=res_file,action="read") + + do i=1,22 + read(15,*) + end do + + open (16, file=analyse_file,position="append",action="write") + n=0 + do i=1,nb_mat + read(15,*)tab(1:7) + read(14,*)mat_file,eig + + if(tab(4)-eig==0)then + write(psb_out_unit,'("on a gagne !")') + n=n+1 + else + write(psb_out_unit,'("perdu")') + write(16,'(F20.2,F20.2)')tab(2),tab(3) + write(16,'(F20.7,F20.7)')tab(4),eig + end if + end do + write(psb_out_unit,'(i20)')n + close(14) + close(15) + close(16) + +end program analyse diff --git a/test/eigen/arnoldi_file.f90 b/test/eigen/arnoldi_file.f90 new file mode 100644 index 00000000..7b0d425d --- /dev/null +++ b/test/eigen/arnoldi_file.f90 @@ -0,0 +1,344 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program arnoldi_file + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_zspmat_type) :: a, aux_a + + ! dense matrices + complex(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + complex(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:) + integer, allocatable :: indexes(:) + type(psb_z_vect_type), allocatable, target :: V(:) + complex(psb_dpk_), pointer :: b_col_glob(:) + type(psb_z_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero, dim_H, alloc_stat + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + real (psb_dpk_) :: norm + complex (psb_dpk_) :: dotprod + integer(psb_ipk_) :: nrhs, nrow, n_row, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='arnoldi_file' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) dim_H + write (psb_out_unit, '("Hessenberg matrix dim is ",i20)') dim_H + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,dim_H) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeig we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + write(psb_out_unit,'("Generating an rhs...")') + write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + write(psb_out_unit,'("Partition type: graph")') + write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + allocate(H(dim_H,dim_H),stat = alloc_stat) + do i=1,dim_h + do j=1,dim_H + H(i,j)=zzero + enddo + enddo + allocate(V(dim_H+1),stat = alloc_stat) + do i=1,dim_H+1 + call psb_geall(V(i),desc_a,info) + call psb_geasb(V(i),desc_a,info) + enddo + + call V(1)%set(zone) + t2 = psb_wtime() - t1 + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + write(psb_out_unit,'(" ")') + end if + + + call psb_barrier(ictxt) + t1 = psb_wtime() + + norm = psb_genrm2(V(1),desc_a,info) + if (iam==psb_root_) write(psb_out_unit,'("norma iniziale :"F20.3)')norm +! H(2,1)=dcmplx(norm,0.0) + H(2,1)=norm*zone + norm = 1/norm + !normalisation of V(1) + call psb_geaxpby(zzero,V(1),norm*zone,V(1),desc_a, info) + do i=2,dim_H+1 + call psb_spmm(zone,a,V(i-1),zzero,V(i),desc_a,info,'n') !we do V(i)=a*V(i-1) + ! Gram-Schmitt's reorthogonalisation + do j=1,i-1 + dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j)) + call psb_geaxpby(-dotprod,V(j),zone,V(i),desc_a, info) !V(i)=V(i)-V(j)*dotprod + if(iam==psb_root_) write(*,'("dotprod : "i5, g20.4,g20.4)') i,real(dotprod),aimag(dotprod) + H(j,i-1)=dotprod + end do + norm = psb_genrm2(V(i),desc_a,info) + if (iam==psb_root_) then + write(psb_out_unit,'("norma finale :"i20,F20.3)')i,norm + write(psb_out_unit,'("")') + end if + if (i .ne. dim_H+1) then + H(i,i-1)=cmplx(norm,0.0) + !H(i,i-1)=norm*zone + endif + norm=1/norm + call psb_geaxpby(zzero,V(i),norm*zone,V(i),desc_a, info) + enddo + do i=2,dim_H + if (iam==psb_root_) write(*,'("basse diagonale de H : "i5, g20.4,g20.4)') i,real(H(i,i-1)),aimag(H(i,i-1)) + enddo + + write(psb_out_unit,'("")') + + allocate(eig(dim_H),work(dim_h),stat = info) + allocate(Z(dim_H,dim_H),stat = alloc_stat) + call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,3*dim_H,info) + + !sort H's eigenvalues + allocate(indexes(1:dim_H)) + call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_) + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + write(psb_out_unit,'("Matrix: ",a)') mtrx_file + write(psb_out_unit,'("Test on : ",i20," processors")') np + write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + + write(*,'("valeurs propres de H : ")') + OPEN(unit=10, file=mtrx_file(1:3)//'.txt') + do i=dim_H/3,dim_H + write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i)) + write(10,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i)) + enddo + CLOSE(10) + + tflops = tflops / (tt2) + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ & + & annz*(psb_sizeof_sp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + bdwdth = times*nbytes/(tt2*1.d6) + + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + do i=1,dim_H + call psb_gefree(V(i),desc_a,info) + end do + DEALLOCATE (H) + DEALLOCATE (eig,V) + DEALLOCATE (work) + DEALLOCATE (Z) +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +end program arnoldi_file + + + + diff --git a/test/eigen/arnoldi_file_real.f90 b/test/eigen/arnoldi_file_real.f90 new file mode 100644 index 00000000..a4befc1e --- /dev/null +++ b/test/eigen/arnoldi_file_real.f90 @@ -0,0 +1,409 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program arnoldi_file + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_dspmat_type) :: a, b, aux_a + + ! dense matrices + real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:) + integer, allocatable :: indexes(:) + type(psb_d_vect_type), allocatable, target :: V(:) + real(psb_dpk_), pointer :: b_col_glob(:) + type(psb_d_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero, dim_H, alloc_stat + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + real (psb_dpk_) :: norm + real (psb_dpk_) :: dotprod + integer(psb_ipk_) :: nrhs, nrow, n_row, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='arnoldi_file' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + !write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + !write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) dim_H + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,dim_H) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeig we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case ('AD') + call adj_read(aux_a,mtrx_file,iunit,desc_a,info) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + ! write(psb_out_unit,'("Generating an rhs...")') + ! write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + !write(psb_out_unit,'("Partition type: graph")') + !write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + !if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call lapl(a,b) + + allocate(H(dim_H,dim_H),stat = alloc_stat) + do i=1,dim_h + do j=1,dim_H + H(i,j)=zzero + enddo + enddo + allocate(V(dim_H+1),stat = alloc_stat) + do i=1,dim_H+1 + call psb_geall(V(i),desc_a,info) + call psb_geasb(V(i),desc_a,info) + enddo + call V(1)%set(done) + t2 = psb_wtime() - t1 + call psb_amx(ictxt, t2) + + ! if (iam==psb_root_) then + !write(psb_out_unit,'(" ")') + !write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + !write(psb_out_unit,'(" ")') + !end if + + call psb_barrier(ictxt) + + norm = psb_norm2(V(1),desc_a,info) + H(2,1)=cmplx(norm,0.0) + norm = 1/norm + !normalisation of V(1) + call psb_geaxpby(dzero,V(1),norm,V(1),desc_a, info) + do i=2,dim_H+1 + call psb_spmm(done,b,V(i-1),dzero,V(i),desc_a,info,'n') !we do V(i)=b*V(i-1) + ! Gram-Schmitt's reorthogonalisation + do j=1,i-1 + dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j)) + call psb_geaxpby(-dotprod,V(j),done,V(i),desc_a, info) !V(i)=V(i)-V(j)*dotprod + H(j,i-1)=cmplx(dotprod,0.0) + end do + norm = psb_norm2(V(i),desc_a,info) + ! write(psb_out_unit,'("norma finale :"i20,F20.3)')i,norm + ! write(psb_out_unit,'("")') + if (i .ne. dim_H+1) then + H(i,i-1)=cmplx(norm,0.0) + endif + norm=1/norm + call psb_geaxpby(dzero,V(i),norm,V(i),desc_a, info) + enddo + + write(psb_out_unit,'("")') + + if (iam==psb_root_) then + allocate(eig(dim_H),work(dim_h),Z(dim_H,dim_H),stat = info) + call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,dim_H,info) + + !sort H's eigenvalues + allocate(indexes(1:dim_H)) + call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_) + end if + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + ! write(psb_out_unit,'("Matrix: ",a)') mtrx_file + ! write(psb_out_unit,'("Test on : ",i20," processors")') np + ! write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + ! write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + + open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE") + write (15, '(a,F20.6,F20.6,F20.4)')mtrx_file,real(eig(dim_H)),real(eig(dim_H-1)),t2 + close(15) + DEALLOCATE (work,eig,Z) + + !write(*,'("valeurs propres de H : ")') + !do i=dim_H/3,dim_H + ! write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i)) + !enddo + + tflops = tflops / (tt2) + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ & + & annz*(psb_sizeof_sp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + bdwdth = times*nbytes/(tt2*1.d6) + + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + do i=1,dim_H + call psb_gefree(V(i), desc_a,info) + enddo + DEALLOCATE (H) + DEALLOCATE (V) +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +contains + + subroutine adj_read (a,filename,iunit,desc_a,info) + type(psb_dspmat_type), intent (inout) :: a + character(len=40) :: filename + integer (psb_ipk_) :: iunit + type(psb_desc_type):: desc_a + integer (psb_ipk_) :: info + + integer(psb_ipk_) :: i,nnzero,nrows + integer (psb_ipk_) :: iError + type(psb_d_coo_sparse_mat) :: acoo + + open(iunit, FILE=filename, STATUS="OLD", ACTION="READ") + read(iunit, *) nrows , nnzero + + call acoo%allocate(nrows,nrows,nnzero) + do i = 1,nnzero + read(iunit, *) acoo%ia(i),acoo%ja(i) + acoo%ia(i)=acoo%ia(i)+1 + acoo%ja(i)=acoo%ja(i)+1 + acoo%val(i)=1.0 + end do + close(UNIT=iunit) + !call psb_spall(a,desc_a,info,nnzero) + !call psb_spins(nnzero, ia, ja, val, a, desc_a, info) + call acoo%set_nzeros(nnzero) + call acoo%fix(info) + call a%mv_from(acoo) + call a%cscnv(info,type='csr') + end subroutine adj_read + + + subroutine lapl(a,b) + type(psb_dspmat_type),intent(in)::a + type(psb_dspmat_type),intent(out)::b + + type(psb_d_coo_sparse_mat) :: acoo + integer(psb_ipk_) :: nz,n,info,i + real(psb_dpk_), allocatable :: K(:) + + call a%cp_to(acoo) + nz=acoo%get_nzeros() + n=a%get_nrows() + allocate(K(n)) + do i=1,n + K(i)=0 + enddo + do i=1,nz + K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i) + acoo%val(i)=-acoo%val(i) + enddo + call acoo%reallocate(nz+n) + do i=1,n + acoo%val(nz+i)=K(i) + acoo%ia(nz+i)= i + acoo%ja(nz+i)= i + enddo + call acoo%set_nzeros(nz+n) + call acoo%fix(info) + call b%mv_from(acoo) + !if(iam==psb_root_) then + ! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")') b%get_nzeros() + !else + ! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros() + !end if + call b%cscnv(info,'CSR') + deallocate (K) + end subroutine lapl + +end program arnoldi_file + + + + diff --git a/test/eigen/d_file_spmv.f90 b/test/eigen/d_file_spmv.f90 new file mode 100644 index 00000000..4db2d0a6 --- /dev/null +++ b/test/eigen/d_file_spmv.f90 @@ -0,0 +1,299 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program d_file_spmv + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_dspmat_type) :: a, aux_a + + ! dense matrices + real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + real(psb_dpk_), pointer :: b_col_glob(:) + type(psb_d_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_), parameter :: times=10 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='d_file_spmv' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeing we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + write(psb_out_unit,'("Generating an rhs...")') + write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + write(psb_out_unit,'("Partition type: graph")') + write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + if (iam==psb_root_) write(psb_out_unit,'("Partition type: default block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call psb_geall(x_col,desc_a,info) + call x_col%set(done) + call psb_geasb(x_col,desc_a,info) + t2 = psb_wtime() - t1 + + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + write(psb_out_unit,'(" ")') + end if + + + call psb_barrier(ictxt) + t1 = psb_wtime() + do i=1,times + call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n') + end do + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + ! FIXME: cache flush needed here + call psb_barrier(ictxt) + tt1 = psb_wtime() + do i=1,times + call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'t') + end do + call psb_barrier(ictxt) + tt2 = psb_wtime() - tt1 + call psb_amx(ictxt,tt2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + write(psb_out_unit,'("Matrix: ",a)') mtrx_file + write(psb_out_unit,'("Test on : ",i20," processors")') np + write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops + flops = flops / (t2) + tflops = tflops / (tt2) + write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2 + write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6 + + write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2 + write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6 + + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_dp + psb_sizeof_int)+& + & annz*(psb_sizeof_dp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + write(psb_out_unit,*) + write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth + bdwdth = times*nbytes/(tt2*1.d6) + write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth + write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt() + + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +end program d_file_spmv + + + + + diff --git a/test/eigen/laplacian.f90 b/test/eigen/laplacian.f90 new file mode 100644 index 00000000..68407f09 --- /dev/null +++ b/test/eigen/laplacian.f90 @@ -0,0 +1,284 @@ +program laplacian + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: mtrx_file, rhs_file + + ! sparse matrices + type(psb_dspmat_type) :: a, b, aux_a + + ! dense matrices + real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + real(psb_dpk_), pointer :: b_col_glob(:) + type(psb_d_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + real(psb_dpk_) :: lambda, lambda2 + real (psb_dpk_) :: norm, precisione + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='laplacian' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + !write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + !write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) precisione + !write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,precisione) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeing we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case ('AD') + call adj_read(aux_a,mtrx_file,iunit,desc_a,info) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + !write(psb_out_unit,'("Generating an rhs...")') + !write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + !if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + !write(psb_out_unit,'("Partition type: graph")') + !write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + !if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call lapl(a,b) + + ! call psb_gefree(b_col, desc_a,info) + !call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_spfree(b, desc_a,info) + call psb_cdfree(desc_a,info) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + + +contains + + subroutine adj_read (a,filename,iunit,desc_a,info) + type(psb_dspmat_type), intent (inout) :: a + character(len=40) :: filename + integer (psb_ipk_) :: iunit + type(psb_desc_type):: desc_a + integer (psb_ipk_) :: info + + integer(psb_ipk_) :: i,nnzero,nrows + integer (psb_ipk_) :: iError + type(psb_d_coo_sparse_mat) :: acoo + + open(iunit, FILE=filename, STATUS="OLD", ACTION="READ") + read(iunit, *) nrows , nnzero + + call acoo%allocate(nrows,nrows,nnzero) + do i = 1,nnzero + read(iunit, *) acoo%ia(i),acoo%ja(i) + acoo%ia(i)=acoo%ia(i)+1 + acoo%ja(i)=acoo%ja(i)+1 + acoo%val(i)=1.0 + end do + close(UNIT=iunit) + call acoo%set_nzeros(nnzero) + call acoo%fix(info) + call a%mv_from(acoo) + call a%cscnv(info,type='csr') + end subroutine adj_read + + + subroutine lapl(a,b) + type(psb_dspmat_type),intent(in)::a + type(psb_dspmat_type),intent(out)::b + + type(psb_d_coo_sparse_mat) :: acoo + integer(psb_ipk_) :: nz,n,info,i + real(psb_dpk_), allocatable :: K(:) + + call a%cp_to(acoo) + nz=acoo%get_nzeros() + n=a%get_nrows() + allocate(K(n)) + do i=1,n + K(i)=0 + enddo + do i=1,nz + K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i) + acoo%val(i)=-acoo%val(i) + enddo + call acoo%reallocate(nz+n) + call acoo%set_dupl(psb_dupl_add_) + do i=1,n + acoo%val(nz+i)=K(i) + acoo%ia(nz+i)= i + acoo%ja(nz+i)= i + enddo + call acoo%set_nzeros(nz+n) + call acoo%fix(info) + do i=1,nz + if(acoo%ja(i)==acoo%ia(i)) then + write(psb_out_unit,'(i10,i10,g20.4)')acoo%ia(i),acoo%ja(i), acoo%val(i) + end if + enddo + !if(iam==psb_root_) then + ! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")')b%get_nzeros() + ! else + ! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros() + !end if + call b%mv_from(acoo) + ! if(iam==psb_root_) then + ! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")')b%get_nzeros() + ! else + ! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n + ! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros() + ! end if + + call b%cscnv(info,'CSR') + deallocate (K) + end subroutine lapl + +end program laplacian diff --git a/test/eigen/pdgenspmv.f90 b/test/eigen/pdgenspmv.f90 new file mode 100644 index 00000000..7df36d37 --- /dev/null +++ b/test/eigen/pdgenspmv.f90 @@ -0,0 +1,308 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +! File: ppde.f90 +! +program pdgenspmv + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer(psb_ipk_) :: idim + + ! miscellaneous + real(psb_dpk_), parameter :: one = 1.d0 + real(psb_dpk_) :: t1, t2, tprec, flops, tflops, tt1, tt2, bdwdth + + ! sparse matrix and preconditioner + type(psb_dspmat_type) :: a + ! descriptor + type(psb_desc_type) :: desc_a + ! dense matrices + type(psb_d_vect_type) :: xv,bv, vtst + real(psb_dpk_), allocatable :: tst(:) + ! blacs parameters + integer(psb_ipk_) :: ictxt, iam, np + + ! solver parameters + integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, precsize, descsize, d2size, annz, nbytes + real(psb_dpk_) :: err, eps + integer(psb_ipk_), parameter :: times=10 + + ! other variables + integer(psb_ipk_) :: info, i + character(len=20) :: name,ch_err + character(len=40) :: fname + + info=psb_success_ + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + if(psb_get_errstatus() /= 0) goto 9999 + name='pde90' + call psb_set_errverbosity(itwo) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + end if + ! + ! get parameters + ! + call get_parms(ictxt,afmt,idim) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& + & a1,a2,a3,b1,b2,b3,c,g,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_gen_pde3d' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 + if (iam == psb_root_) write(psb_out_unit,'(" ")') + + call xv%set(done) + + call psb_barrier(ictxt) + t1 = psb_wtime() + + do i=1,times + call psb_spmm(done,a,xv,dzero,bv,desc_a,info,'n') + end do + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + ! FIXME: cache flush needed here + call psb_barrier(ictxt) + tt1 = psb_wtime() + do i=1,times + call psb_spmm(done,a,xv,dzero,bv,desc_a,info,'t') + end do + call psb_barrier(ictxt) + tt2 = psb_wtime() - tt1 + call psb_amx(ictxt,tt2) + + call psb_amx(ictxt,t2) + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = a%sizeof() + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam == psb_root_) then + flops = 2.d0*times*annz + tflops=flops + write(psb_out_unit,'("Matrix: ell1 ",i0)') idim + write(psb_out_unit,'("Test on : ",i20," processors")') np + write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops + flops = flops / (t2) + tflops = tflops / (tt2) + write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2 + write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6 + + write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2 + write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6 + + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_dp + psb_sizeof_int)+& + & annz*(psb_sizeof_dp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + write(psb_out_unit,*) + write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth + bdwdth = times*nbytes/(tt2*1.d6) + write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth + write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt() + write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize + + end if + + + ! + ! cleanup storage and exit + ! + call psb_gefree(bv,desc_a,info) + call psb_gefree(xv,desc_a,info) + call psb_spfree(a,desc_a,info) + call psb_cdfree(desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='free routine' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + +9999 continue + if(info /= psb_success_) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +contains + ! + ! get iteration parameters from standard input + ! + subroutine get_parms(ictxt,afmt,idim) + integer(psb_ipk_) :: ictxt + character(len=*) :: afmt + integer(psb_ipk_) :: idim + integer(psb_ipk_) :: np, iam + integer(psb_ipk_) :: intbuf(10), ip + + call psb_info(ictxt, iam, np) + + if (iam == 0) then + read(psb_inp_unit,*) afmt + read(psb_inp_unit,*) idim + endif + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,idim) + + if (iam == 0) then + write(psb_out_unit,'("Testing matrix : ell1")') + write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim + write(psb_out_unit,'("Number of processors : ",i0)')np + write(psb_out_unit,'("Data distribution : BLOCK")') + write(psb_out_unit,'(" ")') + end if + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer(psb_ipk_) :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: pde90 methd prec dim & + &[istop itmax itrace]' + write(iout,*)' where:' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' prec : bjac diag none' + write(iout,*)' dim number of points along each axis' + write(iout,*)' the size of the resulting linear ' + write(iout,*)' system is dim**3' + write(iout,*)' istop stopping criterion 1, 2 ' + write(iout,*)' itmax maximum number of iterations [500] ' + write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' >= 1 do tracing every itrace' + write(iout,*)' iterations ' + end subroutine pr_usage + + ! + ! functions parametrizing the differential equation + ! + function b1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b1 + real(psb_dpk_), intent(in) :: x,y,z + b1=1.d0/sqrt(3.d0) + end function b1 + function b2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b2 + real(psb_dpk_), intent(in) :: x,y,z + b2=1.d0/sqrt(3.d0) + end function b2 + function b3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: b3 + real(psb_dpk_), intent(in) :: x,y,z + b3=1.d0/sqrt(3.d0) + end function b3 + function c(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: c + real(psb_dpk_), intent(in) :: x,y,z + c=0.d0 + end function c + function a1(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a1 + real(psb_dpk_), intent(in) :: x,y,z + a1=1.d0/80 + end function a1 + function a2(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a2 + real(psb_dpk_), intent(in) :: x,y,z + a2=1.d0/80 + end function a2 + function a3(x,y,z) + use psb_base_mod, only : psb_dpk_ + real(psb_dpk_) :: a3 + real(psb_dpk_), intent(in) :: x,y,z + a3=1.d0/80 + end function a3 + function g(x,y,z) + use psb_base_mod, only : psb_dpk_, done + real(psb_dpk_) :: g + real(psb_dpk_), intent(in) :: x,y,z + g = dzero + if (x == done) then + g = done + else if (x == dzero) then + g = exp(y**2-z**2) + end if + end function g +end program pdgenspmv diff --git a/test/eigen/power_file.f90 b/test/eigen/power_file.f90 new file mode 100644 index 00000000..64118118 --- /dev/null +++ b/test/eigen/power_file.f90 @@ -0,0 +1,300 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program power_file + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_zspmat_type) :: a, aux_a + + ! dense matrices + complex(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + complex(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + complex(psb_dpk_), pointer :: b_col_glob(:) + type(psb_z_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + complex(psb_dpk_) :: lambda, lambda2 + real (psb_dpk_) :: norm, precisione + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='power_file' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) precisione + write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,precisione) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeing we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + write(psb_out_unit,'("Generating an rhs...")') + write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + write(psb_out_unit,'("Partition type: graph")') + write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call psb_geall(x_col,desc_a,info) + call x_col%set(zone) + call psb_geasb(x_col,desc_a,info) + t2 = psb_wtime() - t1 + + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + write(psb_out_unit,'(" ")') + end if + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_spmm(zone,a,x_col,zzero,b_col,desc_a,info,'n') !we do b_col=a*x_col + norm = psb_norm2(b_col,desc_a,info) + norm = 1/norm + !normalisation of b_col in x_col + call psb_geaxpby(zzero,x_col,norm*zone,b_col,desc_a, info) + x_col=b_col + lambda=0 + lambda2=1000 + do while (abs(lambda-lambda2)>precisione) + !write (psb_out_unit, '("times",i20)')times + times=times+1 + lambda2=lambda + call psb_spmm(zone,a,x_col,zzero,b_col,desc_a,info,'n') + norm = psb_norm2(b_col,desc_a,info) + lambda = psb_gedot(b_col,x_col,desc_a,info) ! lambda = (A*x_col dot x_col) + norm = 1/norm + call psb_geaxpby(zzero,x_col,norm*zone,b_col,desc_a, info) + x_col=b_col + !write (psb_out_unit, '("abs(lambda-lambda2)",F30.20)') abs(lambda-lambda2) + !write (psb_out_unit, '("lambda",F30.20,F30.20)') lambda + ! write (psb_out_unit, '("precisione",F30.20,i20)') precisione, times + end do + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + write(psb_out_unit,'("Matrix: ",a)') mtrx_file + write(psb_out_unit,'("Test on : ",i20," processors")') np + write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops + flops = flops / (t2) + tflops = tflops / (tt2) + write (psb_out_unit, '("The biggest eigenvalue of A is : ",F20.3,F20.3)')lambda + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ & + & annz*(psb_sizeof_sp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + bdwdth = times*nbytes/(tt2*1.d6) + + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +end program power_file + + + + diff --git a/test/eigen/power_file_real.f90 b/test/eigen/power_file_real.f90 new file mode 100644 index 00000000..20f7c5dd --- /dev/null +++ b/test/eigen/power_file_real.f90 @@ -0,0 +1,343 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program power_file_real + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=100) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_dspmat_type) :: a, aux_a + + ! dense matrices + real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + real(psb_dpk_), pointer :: b_col_glob(:) + type(psb_d_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + real(psb_dpk_) :: lambda, lambda2 + real (psb_dpk_) :: norm, precisione + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='power_file' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + !write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + !write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) precisione + !write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,precisione) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeing we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case ('AD') + call adj_read(aux_a,mtrx_file,iunit,desc_a,info) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + !write(psb_out_unit,'("Generating an rhs...")') + !write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + !if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + !write(psb_out_unit,'("Partition type: graph")') + !write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + !if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call psb_geall(x_col,desc_a,info) + call x_col%set(done) + call psb_geasb(x_col,desc_a,info) + t2 = psb_wtime() - t1 + + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + !write(psb_out_unit,'(" ")') + !write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + !write(psb_out_unit,'(" ")') + end if + + call psb_barrier(ictxt) + call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n') !we do b_col=a*x_col + norm = psb_norm2(b_col,desc_a,info) + norm = 1/norm + !normalisation of b_col in x_col + call psb_geaxpby(dzero,x_col,norm,b_col,desc_a, info) + x_col=b_col + lambda=0 + lambda2=1000 + do while (abs(lambda-lambda2)>precisione .AND. times<300) + !if(iam==psb_root_) write (psb_out_unit,'(F20.6,i20)')abs(lambda-lambda2), times + times=times+1 + lambda2=lambda + call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n') + norm = psb_norm2(b_col,desc_a,info) + lambda = psb_gedot(b_col,x_col,desc_a,info) ! lambda = (A*x_col dot x_col) + norm = 1/norm + call psb_geaxpby(dzero,x_col,norm,b_col,desc_a, info) + x_col=b_col + end do + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + ! write(psb_out_unit,'("Matrix: ",a)') mtrx_file + ! write(psb_out_unit,'("Test on : ",i20," processors")') np + ! write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + ! write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + !write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + ! write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops + ! write(psb_out_unit,'("Time for ",i0," iteration (s) : ",F20.3)')times, t2 + ! write(psb_out_unit,'("Time per iteration (ms) : ",F20.3)') t2*1.d3/(1.d0*times) + open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE") + write (15, '(a,F20.6,F20.4)')mtrx_file,lambda,t2 + ! write (psb_out_unit, '("entropy vect propre :",g20.4)')entropy(x_col) + close(15) + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +contains + subroutine adj_read (a,filename,iunit,desc_a,info) + type(psb_dspmat_type), intent (inout) :: a + character(len=100) :: filename + integer (psb_ipk_) :: iunit + type(psb_desc_type):: desc_a + integer (psb_ipk_) :: info + + integer(psb_ipk_) :: i,nnzero,nrows + integer (psb_ipk_) :: iError + !integer(psb_ipk_), allocatable :: ia(:),ja(:) + !real (psb_dpk_), allocatable :: val(:) + type(psb_d_coo_sparse_mat) :: acoo + + open(iunit, FILE=filename, STATUS="OLD", ACTION="READ") + read(iunit, *) nrows , nnzero + + call acoo%allocate(nrows,nrows,nnzero) + do i = 1,nnzero + read(iunit, *) acoo%ia(i),acoo%ja(i) + acoo%ia(i)=acoo%ia(i)+1 + acoo%ja(i)=acoo%ja(i)+1 + acoo%val(i)=1.0 + end do + close(UNIT=iunit) + !call psb_spall(a,desc_a,info,nnzero) + !call psb_spins(nnzero, ia, ja, val, a, desc_a, info) + call acoo%set_nzeros(nnzero) + call acoo%fix(info) + call a%mv_from(acoo) + call a%cscnv(info,type='csr') + end subroutine adj_read + + + real(psb_dpk_) function entropy(vect) + type(psb_d_vect_type), intent(inout)::vect + + real(psb_dpk_), allocatable :: v(:) + integer(psb_ipk_) :: nrows + + nrows=vect%get_nrows() + allocate(v(nrows)) + v=vect%get_vect() + entropy=0 + do i=1,nrows + entropy=entropy+v(i)*v(i)*log(v(i)*v(i)) + end do + return + end function entropy +end program power_file_real + + + + diff --git a/test/eigen/s_file_spmv.f90 b/test/eigen/s_file_spmv.f90 new file mode 100644 index 00000000..fddff2e5 --- /dev/null +++ b/test/eigen/s_file_spmv.f90 @@ -0,0 +1,298 @@ +!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program s_file_spmv + use psb_base_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices + type(psb_sspmat_type) :: a, aux_a + + ! dense matrices + real(psb_spk_), allocatable, target :: aux_b(:,:), d(:) + real(psb_spk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + real(psb_spk_), pointer :: b_col_glob(:) + type(psb_s_vect_type) :: b_col, x_col, r_col + + + ! communications data structure + type(psb_desc_type):: desc_a + + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr + integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes + real(psb_spk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_), parameter :: times=10 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='s_file_spmv' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeing we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + write(psb_out_unit,'("Generating an rhs...")') + write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + write(psb_out_unit,'("Partition type: graph")') + write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call psb_geall(x_col,desc_a,info) + call x_col%set(sone) + call psb_geasb(x_col,desc_a,info) + t2 = psb_wtime() - t1 + + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + write(psb_out_unit,'(" ")') + end if + + + call psb_barrier(ictxt) + t1 = psb_wtime() + do i=1,times + call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'n') + end do + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + ! FIXME: cache flush needed here + call psb_barrier(ictxt) + tt1 = psb_wtime() + do i=1,times + call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'t') + end do + call psb_barrier(ictxt) + tt2 = psb_wtime() - tt1 + call psb_amx(ictxt,tt2) + + nr = desc_a%get_global_rows() + annz = a%get_nzeros() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + write(psb_out_unit,'("Matrix: ",a)') mtrx_file + write(psb_out_unit,'("Test on : ",i20," processors")') np + write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops + flops = flops / (t2) + tflops = tflops / (tt2) + write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2 + write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6 + + write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2 + write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times) + write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6 + + ! + ! This computation is valid for CSR + ! + nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ & + & annz*(psb_sizeof_sp + psb_sizeof_int) + bdwdth = times*nbytes/(t2*1.d6) + write(psb_out_unit,*) + write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth + bdwdth = times*nbytes/(tt2*1.d6) + write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth + + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_cdfree(desc_a,info) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +end program s_file_spmv + + + + + diff --git a/test/eigen/shift_invert_real.f90 b/test/eigen/shift_invert_real.f90 new file mode 100644 index 00000000..a12befcf --- /dev/null +++ b/test/eigen/shift_invert_real.f90 @@ -0,0 +1,471 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +program shift_invert + use psb_base_mod + use psb_prec_mod + use psb_krylov_mod + use psb_util_mod + implicit none + + ! input parameters + character(len=40) :: kmethd, ptype, mtrx_file, rhs_file + + ! sparse matrices and preconditioner + type(psb_dspmat_type) :: a, aux_a, id, a2, b + type(psb_d_csr_sparse_mat)::acsr, icsr + type(psb_dprec_type) :: prec + + ! dense matrices + real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) + complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:) + + real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:) + real(psb_dpk_), pointer :: b_col_glob(:) + type(psb_d_vect_type) :: b_col, x_col, r_col + type (psb_d_vect_type), allocatable, target :: V(:) + integer(psb_ipk_), allocatable , save :: indexes(:) + + + ! communications data structure + type(psb_desc_type):: desc_a + integer(psb_ipk_) :: ictxt, iam, np + + ! solver paramters + integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,& + & methd, istopc, irst, nr, sort, annz + integer(psb_long_int_k_) :: amatsize, descsize, nbytes + real(psb_dpk_) :: err, eps,cond + + character(len=5) :: afmt + character(len=20) :: name, ch_err + character(len=2) :: filefmt + integer(psb_ipk_), parameter :: iunit=12 + integer(psb_ipk_) :: times=0 + integer(psb_ipk_) :: iparm(20) + + ! other variables + integer(psb_ipk_) :: i,info,j,m_problem + integer(psb_ipk_) :: internal, m,ii,nnzero,ji + real(psb_dpk_) :: t1, t2, r_amax, b_amax,& + &scale,resmx,resmxp, flops, bdwdth + real(psb_dpk_) :: tt1, tt2, tflops, tprec, sigma,dotprod,norm + integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne,dim_H + integer(psb_ipk_), allocatable :: ivg(:), ipv(:) + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + + + name='shift_invert_real' + if(psb_get_errstatus() /= 0) goto 9999 + info=psb_success_ + call psb_set_errverbosity(2) + ! + ! Hello world + ! + if (iam == psb_root_) then + write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ + write(*,*) 'This is the ',trim(name),' sample program' + read(psb_inp_unit,*) mtrx_file + read(psb_inp_unit,*) filefmt + read(psb_inp_unit,*) ipart + read(psb_inp_unit,*) ptype + read(psb_inp_unit,*) kmethd + read(psb_inp_unit,*) dim_H + read(psb_inp_unit,*) sigma + read(psb_inp_unit,*) eps + end if + call psb_bcast(ictxt,mtrx_file) + call psb_bcast(ictxt,filefmt) + call psb_bcast(ictxt,ipart) + call psb_bcast(ictxt,dim_H) + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,ptype) + call psb_bcast(ictxt,sigma) + call psb_bcast(ictxt,eps) + rhs_file = 'NONE' + afmt = 'CSR' + call psb_barrier(ictxt) + t1 = psb_wtime() + ! read the input matrix to be processed and (possibly) the rhs + nrhs = 1 + + if (iam==psb_root_) then + select case(psb_toupper(filefmt)) + case('MM') + ! For Matrix Market we have an input file for the matrix + ! and an (optional) second file for the RHS. + call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file) + if (info == psb_success_) then + if (rhs_file /= 'NONE') then + call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file) + end if + end if + + case ('HB') + ! For Harwell-Boeig we have a single file which may or may not + ! contain an RHS. + call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file) + + case ('AD') + call adj_read(aux_a,mtrx_file,iunit,desc_a,info) + + case default + info = -1 + write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt + end select + if (info /= psb_success_) then + write(psb_err_unit,*) 'Error while reading input matrix ' + call psb_abort(ictxt) + end if + + m_problem = aux_a%get_nrows() + annz=aux_a%get_nzeros() + call psb_bcast(ictxt,m_problem) + + ! At this point aux_b may still be unallocated + if (psb_size(aux_b,dim=1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(psb_err_unit,'("Ok, got an rhs ")') + b_col_glob =>aux_b(:,1) + else + ! write(psb_out_unit,'("Generating an rhs...")') + ! write(psb_out_unit,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + + b_col_glob => aux_b(:,1) + do i=1, m_problem + b_col_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + else + + call psb_bcast(ictxt,m_problem) + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + endif + b_col_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_col_glob(1:m_problem)) + + end if + + ! switch over different partition types + if (ipart == 0) then + call psb_barrier(ictxt) + if (iam==psb_root_) write(psb_out_unit,'("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 psb_matdist(aux_a, a,ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else if (ipart == 2) then + if (iam==psb_root_) then + !write(psb_out_unit,'("Partition type: graph")') + !write(psb_out_unit,'(" ")') + ! write(psb_err_unit,'("Build type: graph")') + call build_mtpart(aux_a,np) + + endif + call psb_barrier(ictxt) + call distr_mtpart(psb_root_,ictxt) + call getv_mtpart(ivg) + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg) + + else + if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")') + call psb_matdist(aux_a, a, ictxt, & + & desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block) + end if + + call lapl(a,b) + + ! + !a2=a-sigma*id + ! + !sigma = 100 + call csshift(b,a2,sigma,ictxt) + + ! + ! prepare the preconditioner. + ! + !if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ", a) ') ptype + call psb_precinit(prec,ptype,info) + + call psb_barrier(ictxt) + !t1 = psb_wtime() + call psb_precbld(a2,desc_a,prec,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + tprec = psb_wtime()-t1 + if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec + + allocate(H(dim_H,dim_H)) + allocate (V(dim_H+1)) + do i=1,dim_H+1 + call psb_geall(V(i),desc_a,info) + call psb_geasb(V(i),desc_a,info) + enddo + call V(1)%set(done) + call psb_amx(ictxt, t2) + + !if (iam==psb_root_) then + ! write(psb_out_unit,'(" ")') + ! write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2 + ! write(psb_out_unit,'(" ")') + !end if + + call psb_barrier(ictxt) +t2 = psb_wtime()-t1-tprec + if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')t2 + + + norm = psb_norm2(V(1),desc_a, info) + H(2,1)=cmplx(norm,0.0) + norm = 1/norm + call psb_geaxpby(dzero,V(1),norm,V(1),desc_a,info) + do i=2,dim_H+1 + !A*V(i)=V(i-1) + call psb_krylov(kmethd,a2,prec,V(i-1),V(i),eps,desc_a,info,& + & itmax=1000,iter=iter,err=err,itrace=-1,istop=2,irst=002) + !if (iam==psb_root_) write(*,'("iter : "i20)') iter + ! Gram-Schmitt's reorthogonalisation + do j=1,i-1 + dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j)) + call psb_geaxpby(-dotprod,V(j),done,V(i),desc_a, info)!V(i)=V(i)-V(j)*dotprod + H(j,i-1)=cmplx(dotprod,0.0) + end do + norm = psb_norm2(V(i),desc_a,info) + if (i .ne. dim_H+1) then + H(i,i-1)=cmplx(norm,0.0) + endif + norm=1/norm + call psb_geaxpby(dzero,V(i),norm,V(i),desc_a, info) + enddo + + t2=psb_wtime()-t1-tprec + if (iam==psb_root_) write (psb_out_unit,'("temps de arnoldi : " ,es12.5 )') t2 + if(iam==psb_root_) then + allocate(eig(dim_H),work(dim_h),Z(dim_H,dim_H),stat = info) + call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,dim_H,info) + + !sort H's eigenvalues + allocate(indexes(1:dim_H)) + call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_) + do i=1,dim_H + eig (i) = cmplx(sigma,0.0)+1/eig(i) + !write(psb_out_unit, '("eig(i), i", g20.4, i10)')real(eig(i)),i + enddo + end if + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + nr = desc_a%get_global_rows() + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + call psb_sum(ictxt,annz) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + + if (iam==psb_root_) then + flops = 2.d0*times*annz + tflops=flops + !write(psb_out_unit,'("Matrix: ",a)') mtrx_file + !write(psb_out_unit,'("Test on : ",i20," processors")') np + !write(psb_out_unit,'("Size of matrix : ",i20," ")') nr + !write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz + !write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize + !write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops + + !write(*,'("eigenvalues near from ", g20.4," : ")') sigma + !do i=dim_H/3,dim_H + ! write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i)) + !enddo + open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE") + write (15,'(F20.6,F20.6,F20.4)')real(eig(dim_H-1)),real(eig(dim_H)),t2 + close(15) + DEALLOCATE (work,eig,Z) + end if + + call psb_gefree(b_col, desc_a,info) + call psb_gefree(x_col, desc_a,info) + call psb_spfree(a, desc_a,info) + call psb_spfree(a2, desc_a,info) + call psb_spfree(b, desc_a,info) + call psb_cdfree(desc_a,info) + do i=1,dim_H + call psb_gefree(V(i), desc_a,info) + enddo + DEALLOCATE (H) + DEALLOCATE (V) + +9999 continue + if(info /= 0) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + + +contains + + subroutine csshift(a,b,sigma,ictx) + type(psb_dspmat_type), intent(in) :: a + type(psb_dspmat_type), intent(out) :: b + real(psb_dpk_) :: sigma + integer(psb_ipk_)::ictx + + type(psb_d_coo_sparse_mat) :: acoo + integer(psb_ipk_) :: nz,n,info,i + + call a%cp_to(acoo) + if (sigma/=0.0) then + nz=acoo%get_nzeros() + n=a%get_nrows() + call acoo%reallocate(nz+n) + call acoo%set_dupl(psb_dupl_add_) + do i=1,n + acoo%val(nz+i)=-sigma + acoo%ia(nz+i)= i + acoo%ja(nz+i)= i + enddo + call acoo%set_nzeros(nz+n) + call acoo%fix(info) + ! do i=1,nz + ! if(acoo%ja(i)==acoo%ia(i)) then + ! write(psb_out_unit,'(i10,i10,g20.4)')acoo%ja(i),i, acoo%val(i) + ! end if + ! enddo + !write(psb_out_unit,'("autant de nzeros apres fix ?",i10)') acoo%get_nzeros()-nz + end if + call b%mv_from(acoo) + call b%cscnv(info,'CSR') + end subroutine csshift + + subroutine adj_read (a,filename,iunit,desc_a,info) + type(psb_dspmat_type), intent (inout) :: a + character(len=40) :: filename + integer (psb_ipk_) :: iunit + type(psb_desc_type):: desc_a + integer (psb_ipk_) :: info + + integer(psb_ipk_) :: i,nnzero,nrows + integer (psb_ipk_) :: iError + type(psb_d_coo_sparse_mat) :: acoo + + open(iunit, FILE=filename, STATUS="OLD", ACTION="READ") + read(iunit, *) nrows , nnzero + + call acoo%allocate(nrows,nrows,nnzero) + do i = 1,nnzero + read(iunit, *) acoo%ia(i),acoo%ja(i) + acoo%ia(i)=acoo%ia(i)+1 + acoo%ja(i)=acoo%ja(i)+1 + acoo%val(i)=1.0 + end do + close(UNIT=iunit) + !call psb_spall(a,desc_a,info,nnzero) + !call psb_spins(nnzero, ia, ja, val, a, desc_a, info) + call acoo%set_nzeros(nnzero) + call acoo%fix(info) + call a%mv_from(acoo) + call a%cscnv(info,type='csr') + end subroutine adj_read + +subroutine lapl(a,b) + type(psb_dspmat_type),intent(in)::a + type(psb_dspmat_type),intent(out)::b + + type(psb_d_coo_sparse_mat) :: acoo + integer(psb_ipk_) :: nz,n,info,i + real(psb_dpk_), allocatable :: K(:) + + call a%cp_to(acoo) + nz=acoo%get_nzeros() + n=a%get_nrows() + allocate(K(n)) + do i=1,n + K(i)=0 + enddo + do i=1,nz + K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i) + acoo%val(i)=-acoo%val(i) + enddo + call acoo%reallocate(nz+n) + call acoo%set_dupl(psb_dupl_add_) + do i=1,n + acoo%val(nz+i)=K(i) + acoo%ia(nz+i)= i + acoo%ja(nz+i)= i + enddo + call acoo%set_nzeros(nz+n) + call acoo%fix(info) + do i=1,nz + ! if(acoo%ja(i)==acoo%ia(i)) then + ! write(psb_out_unit,'(i10,i10,g20.4)')acoo%ia(i),acoo%ja(i),acoo%val(i) + ! end if + enddo + call b%mv_from(acoo) + call b%cscnv(info,'CSR') + deallocate (K) + end subroutine lapl + +end program shift_invert + + + + diff --git a/test/eigen/test_chseqr.f90 b/test/eigen/test_chseqr.f90 new file mode 100644 index 00000000..622f911f --- /dev/null +++ b/test/eigen/test_chseqr.f90 @@ -0,0 +1,20 @@ +program test_chseqr + +complex , dimension (1:3,1:3) :: H,Z,Vr +complex , dimension (1:6,1:6) :: Rwork +complex, dimension (1:3) :: eing,work +integer :: info, N +N=3 +H = reshape((/ (1.0,0),(0,0), (2,0), (2,0) , (3,0) , (-4,0) , (0,0) , (0,0) , (2,0) /), shape(H)) +do i=1,N + do j=1,N + write(*,'("H : "i5,i5,2x,g20.4,g20.4)')i,j,real(H(i,j)),aimag(H(i,j)) + end do +end do +call CHSEQR('E','N',N,1,N,H,N,eing,Z,N,work,N,info) +!call CGEEV('N','N',N,H,N,eing,Z,N,Vr,N,work,3*N,Rwork,info) +do i=1,N + write(*,'("valeur propre de H : "i5, i5, g20.4,g20.4)')info, i,real(eing(i)),aimag(eing(i)) +enddo + +end program test_chseqr diff --git a/test/kernel/Makefile b/test/kernel/Makefile index 4d26580f..314f63e5 100644 --- a/test/kernel/Makefile +++ b/test/kernel/Makefile @@ -5,7 +5,7 @@ BASEDIR=../.. INCDIR=$(BASEDIR)/include/ include $(INCDIR)/Make.inc.psblas LIBDIR=$(BASEDIR)/lib/ -PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base +PSBLAS_LIB= -L$(LIBDIR) -leigen -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base LDLIBS=$(PSBLDLIBS) FINCLUDES=$(FMFLAG)$(INCDIR) $(FMFLAG). @@ -47,7 +47,7 @@ arnoldi_file: $(ARNOBJS) $(F90LINK) $(LOPT) $(ARNOBJS) -o arnoldi_file $(PSBLAS_LIB) $(LDLIBS) /bin/mv arnoldi_file $(EXEDIR) -arnoldi_file_real: $(ARNROBJS) +arnoldi_file_real: $(ARNROBJS) . $(F90LINK) $(LOPT) $(ARNROBJS) -o arnoldi_file_real $(PSBLAS_LIB) $(LDLIBS) /bin/mv arnoldi_file_real $(EXEDIR) diff --git a/test/kernel/arnoldi_file.f90 b/test/kernel/arnoldi_file.f90 index 9110bb32..7b0d425d 100644 --- a/test/kernel/arnoldi_file.f90 +++ b/test/kernel/arnoldi_file.f90 @@ -223,6 +223,7 @@ program arnoldi_file call psb_geall(V(i),desc_a,info) call psb_geasb(V(i),desc_a,info) enddo + call V(1)%set(zone) t2 = psb_wtime() - t1 call psb_amx(ictxt, t2) @@ -277,7 +278,7 @@ program arnoldi_file !sort H's eigenvalues allocate(indexes(1:dim_H)) - call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_) + call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_) call psb_barrier(ictxt) t2 = psb_wtime() - t1