From a9722dc049d558e467a89f49ee39cc8e8b5b40a7 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sat, 19 Jul 2008 15:22:51 +0000 Subject: [PATCH] mld2p4: docs/pdf/building.tex docs/pdf/distribution.tex docs/pdf/license.tex docs/pdf/userinterface.tex examples examples/fileread examples/fileread/Makefile examples/fileread/data_input.f90 examples/fileread/mld_dexample_1lev.f90 examples/fileread/mld_dexample_ml.f90 examples/fileread/runs examples/pdegen examples/pdegen/Makefile examples/pdegen/data_input.f90 examples/pdegen/mld_dexample_1lev.f90 examples/pdegen/mld_dexample_ml.f90 examples/pdegen/runs test/fileread/Makefile test/fileread/df_bench.f90 test/fileread/runs/drt.sh test/fileread/zf_bench.f90 test/pargen/Makefile test/pargen/ppde.f90 test/pargen/spde.f90 Added exampels directory. Cleaned up/aligned test/pargen. Moved out Xf_bench to psb_bench. --- docs/pdf/building.tex | 58 ++- docs/pdf/distribution.tex | 13 +- docs/pdf/license.tex | 4 +- docs/pdf/userinterface.tex | 13 +- examples/fileread/Makefile | 40 ++ examples/fileread/data_input.f90 | 91 ++++ examples/fileread/mld_dexample_1lev.f90 | 309 +++++++++++ examples/fileread/mld_dexample_ml.f90 | 348 +++++++++++++ examples/pdegen/Makefile | 40 ++ examples/pdegen/data_input.f90 | 91 ++++ examples/pdegen/mld_dexample_1lev.f90 | 598 ++++++++++++++++++++++ examples/pdegen/mld_dexample_ml.f90 | 636 +++++++++++++++++++++++ test/fileread/Makefile | 14 +- test/fileread/df_bench.f90 | 648 ------------------------ test/fileread/runs/drt.sh | 84 --- test/fileread/zf_bench.f90 | 594 ---------------------- test/pargen/Makefile | 12 +- test/pargen/ppde.f90 | 80 +-- test/pargen/spde.f90 | 91 +--- 19 files changed, 2274 insertions(+), 1490 deletions(-) create mode 100644 examples/fileread/Makefile create mode 100644 examples/fileread/data_input.f90 create mode 100644 examples/fileread/mld_dexample_1lev.f90 create mode 100644 examples/fileread/mld_dexample_ml.f90 create mode 100644 examples/pdegen/Makefile create mode 100644 examples/pdegen/data_input.f90 create mode 100644 examples/pdegen/mld_dexample_1lev.f90 create mode 100644 examples/pdegen/mld_dexample_ml.f90 delete mode 100644 test/fileread/df_bench.f90 delete mode 100755 test/fileread/runs/drt.sh delete mode 100644 test/fileread/zf_bench.f90 diff --git a/docs/pdf/building.tex b/docs/pdf/building.tex index 7fc35840..e3139507 100644 --- a/docs/pdf/building.tex +++ b/docs/pdf/building.tex @@ -1,13 +1,63 @@ \section{Configuring and Building MLD2P4\label{sec:building}} \markboth{\textsc{MLD2P4 User's and Reference Guide}} {\textsc{\ref{sec:building} Configuring and Building MLD2P4}} - - uso di GNU autoconf e automake, con riferimento alla compilazione delle varie versioni - (reale/complessa, singola/doppia precisione) \\ - - software di base (MPI, BLACS, BLAS, PSBLAS - specificare versioni)\\ +To build MLD2P4 it is necessary to set up a Makefile with appropriate +values for your system; this is done by means of the \verb|configure| +script. The distribution also includes the autoconf and automake +sources employed to generate the script, but this should not normally +be needed to build the software. + +MLD2P4 is implemented almost entirely in Fortran~95, with some +interfaces to external libraries in C; we require the Fortran compiler +to support the Fortran~95 standard plus the extension TR15581, which +enhances the usability of \verb|ALLOCATABLE| variables. Most modern +Fortran compilers, including the GNU Fortran compiler, support this +language level. The software defines data types and interfaces for +real and complex data, in both single and double precision. + +\subsection{Prerequisites} +The following base libraries are needed: +\begin{description} +\item[BLAS] The Basic Linear Algebra subprograms. Many vendors provide + optimized versions; if no vendor version is available for a given + platform, the ATLAS software \verb!http://www.netlib.org/atlas! + may be employed. The reference BLAS from Netlib + \verb|htt://www.netlib.org/blas| are meant to define the standard + behaviour of the BLAS interface, so they not optimized for any + particular plaftorm, and they should only be used as a last + resort; note that BLAS computation form a relatively small part of + the MLD2P4/PSBLAS computations, except when using preconditioners + based on the UMFPACK or SuperLU third party libraries. +\item[MPI] A version of MPI is available on most high performance + computing system; we only require version 1.1. +\item[BLACS] The Basic Linear Algebra Communication Subroutines are + available in source form from \verb|http://www.netlib.org/blacs|; + some vendors include them in their parallel computing + support libraries. +\end{description} + +The MLD2P4 software requires PSBLAS version 2.3 (at least), available +from \verb|http://www.ce.uniroma2.it/psblas|; indeed, all the +prerequisites listed so fare are also prerequisites of PSBLAS. Please +note that to build the MLD2P4 library it is necessary to get access to +the source PSBLAS directory used to build the version under use; after +the build process completes, only the compiled form of the library is +necessary to build user applications. + +Please note that all the libraries listed so fare (BLAS, MPI, BLACS, +PSBLAS) must have Fortran interfaces compatible with the MLD2P4; +usually this means that they should all be built with the same +compiler. + + +\subsection{Optional third party libraries} + - software `third party" (UMFPACK, SuperLU, SuperLUdist - specificare versioni e opzioni di configure, dire che UMFPACK e SuperLU sono richiesti, rispettivamente, dalle versioni in singola ed in doppia precisione.)\\ - - sistemi operativi e compilatori su cui MLD2P4 e' stato implementato con successo \\ +\subsection{Configuration options} + - sistemi operativi e compilatori su cui MLD2P4 e' stato + implementato con successo \\ - sono previste opzioni di configurazione per il debugging o per il profiling? \\ - albero delle directory generato al momento dell'installazione, con brevissima descrizione del contenuto delle directory - aggiungere la directory \texttt{examples}\\ diff --git a/docs/pdf/distribution.tex b/docs/pdf/distribution.tex index 8bc01b41..c8911fdc 100644 --- a/docs/pdf/distribution.tex +++ b/docs/pdf/distribution.tex @@ -3,5 +3,14 @@ {\textsc{\ref{sec:distribution} Code Distribution}} \noindent -- sito web sul quale e' disponibile il software\\ -- modalita' di uso con cenno alla licenza (i dettagli sulla licenza vanno in appendice) \ No newline at end of file +MLD2P4 is available from our project web site +\begin{quotation} +\tt http://www.mld2p4.it +\end{quotation} +where you will also find contact points for further information and +bug reports. + +The software is available under a modified BSD license, as specified +in appendix~\ref{sec:license}; please note that some of the optional +third party libraries may be licensed under a different and more +stringent license. diff --git a/docs/pdf/license.tex b/docs/pdf/license.tex index 8f305b24..005ffa08 100644 --- a/docs/pdf/license.tex +++ b/docs/pdf/license.tex @@ -1,6 +1,6 @@ -\section{License\label{sec:distribution}} +\section{License\label{sec:license}} \markboth{\textsc{MLD2P4 User's and Reference Guide}} - {\textsc{\ref{sec:distribution} License}} + {\textsc{\ref{sec:license} License}} The MLD2P4 is freely distributable under the following copyright terms: {\small diff --git a/docs/pdf/userinterface.tex b/docs/pdf/userinterface.tex index 02ed7a26..4813103c 100644 --- a/docs/pdf/userinterface.tex +++ b/docs/pdf/userinterface.tex @@ -57,8 +57,9 @@ according to the preconditioner type chosen by the user. must be chosen according to the real/complex, single/double precision version of MLD2P4 under use.\\ \verb|ptype| & \verb|character(len=*), intent(in)|.\\ - & The type of preconditioner. Its values are specified in Table~\ref{tab:precinit}.\\ - Note that the uppercase and lowercase letters can be used indifferently.\\ + & The type of preconditioner. Its values are specified + in Table~\ref{tab:precinit}.\\ + & Note that the strings are case insensitive.\\ \verb|info| & \verb|integer, intent(out)|.\\ & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\ \verb|nlev| & \verb|integer, optional, intent(in)|.\\ @@ -97,7 +98,7 @@ contained in \verb|val|. values and the corresponding data types is given in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}. When the value is of type \verb|character(len=*)|, - uppercase and lowercase letters can be used indifferently.\\ + it is also treated as case insensitive.\\ \verb|info| & \verb|integer, intent(out)|.\\ & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\ @@ -121,10 +122,10 @@ can be logically divided into four groups, i.e.\ parameters defining \item the aggregation algorithm; \item the coarse-space correction at the coarsest level. \end{enumerate} -A list of the parameters that can be set, along with their allowed and +A list of the parameters that can be set, along with allowed and default values, is given in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}. -For a better understanding of the meaning of the parameters, the user -is referred to Section \ref{sec:background}. +For a detailed description of the meaning of the parameters, please +see Section~\ref{sec:background}. % %Note that the routine allows to set different features of the %preconditioner at each level through the use of \verb|ilev|. diff --git a/examples/fileread/Makefile b/examples/fileread/Makefile new file mode 100644 index 00000000..affefbf8 --- /dev/null +++ b/examples/fileread/Makefile @@ -0,0 +1,40 @@ +MLDDIR=../.. +include $(MLDDIR)/Make.inc +PSBDIR=$(PSBLASDIR)/lib/ +MLDLIBDIR=$(MLDDIR)/lib +MLD_LIB=-L$(MLDLIBDIR) -lmld_krylov -lmld_prec +PSBLAS_LIB= -L$(PSBDIR) -lpsb_util -lpsb_base +FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBDIR) $(FIFLAG). + +DMOBJS=mld_dexample_ml.o data_input.o +D1OBJS=mld_dexample_1lev.o data_input.o + +EXEDIR=./runs + +all: mld_dexample_ml mld_dexample_1lev + +mld_dexample_ml: $(DMOBJS) + $(F90LINK) $(LINKOPT) $(DMOBJS) -o mld_dexample_ml \ + $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv mld_dexample_ml $(EXEDIR) + +mld_dexample_1lev: $(D1OBJS) + $(F90LINK) $(LINKOPT) $(D1OBJS) -o mld_dexample_1lev \ + $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv mld_dexample_1lev $(EXEDIR) + +mld_dexample_ml.o: data_input.o +mld_dexample_1lev.o: data_input.o + +.f90.o: + $(MPF90) $(F90COPT) $(FINCLUDES) -c $< + +clean: + /bin/rm -f $(DMOBJS) $(D1OBJS) \ + *$(.mod) $(EXEDIR)/mld_dexample_ml $(EXEDIR)/mld_dexample_1lev + +lib: + (cd ../../; make library) +verycleanlib: + (cd ../../; make veryclean) + diff --git a/examples/fileread/data_input.f90 b/examples/fileread/data_input.f90 new file mode 100644 index 00000000..32aae9fb --- /dev/null +++ b/examples/fileread/data_input.f90 @@ -0,0 +1,91 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.2) +!!$ +!!$ (C) Copyright 2008 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +module data_input + + interface read_data + module procedure read_char, read_int,& + & read_double, read_single + end interface read_data + +contains + + subroutine read_char(val,file) + character(len=*), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),'(a)') val + end subroutine read_char + subroutine read_int(val,file) + integer, intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_int + subroutine read_single(val,file) + use psb_base_mod + real(psb_spk_), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_single + subroutine read_double(val,file) + use psb_base_mod + real(psb_dpk_), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_double +end module data_input + diff --git a/examples/fileread/mld_dexample_1lev.f90 b/examples/fileread/mld_dexample_1lev.f90 new file mode 100644 index 00000000..a91d5444 --- /dev/null +++ b/examples/fileread/mld_dexample_1lev.f90 @@ -0,0 +1,309 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.2) +!!$ +!!$ (C) Copyright 2008 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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: mld_dexample_ml.f90 +! +! This sample program solves a linear system by using BiCGStab preconditioned by +! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 6.1 +! of the MLD2P4 User's and Reference Guide. +! +! The matrix and the rhs are read from files (if an rhs is not available, the +! unit rhs is set). +! +program mld_dexample_ml + use psb_base_mod + use mld_prec_mod + use psb_krylov_mod + use psb_util_mod + use data_input + + implicit none + +! input parameters + character(len=40) :: mtrx_file, rhs_file + +! sparse matrices + type(psb_dspmat_type) :: A, aux_A + +! descriptor of sparse matrices + type(psb_desc_type):: desc_A + +! preconditioner + type(mld_dprec_type) :: P + +! right-hand side, solution and residual vectors + real(psb_dpk_), allocatable , save :: b(:), x(:), r(:), & + & x_glob(:), r_glob(:) + real(psb_dpk_), allocatable, target :: aux_b(:,:) + real(psb_dpk_), pointer :: b_glob(:) + +! solver and preconditioner parameters + real(psb_dpk_) :: tol, err + integer :: itmax, iter, istop + integer :: nlev + +! parallel environment parameters + integer :: ictxt, iam, np + +! other variables + integer :: i,info,j,m_problem,amatsize,descsize,precsize + integer :: ierr, ircode + real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp + character(len=20) :: name + +! initialize the parallel environment + + 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='mld_dexample_ml' + if(psb_get_errstatus() /= 0) goto 9999 + info=0 + call psb_set_errverbosity(2) + +! get parameters + + call get_parms(ictxt,mtrx_file,rhs_file,itmax,tol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + +! read and assemble the matrix A and the right-hand side b +! using PSBLAS routines for sparse matrix / vector management + + if (iam==psb_root_) then + call read_mat(mtrx_file, aux_A, ictxt) + + m_problem = aux_A%m + call psb_bcast(ictxt,m_problem) + + if(rhs_file /= 'NONE') then + ! reading an rhs + call read_rhs(rhs_file,aux_b,ictxt) + end if + + if (psb_size(aux_b,1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(0,'("Ok, got an rhs ")') + b_glob =>aux_b(:,1) + else + write(*,'("Generating an rhs...")') + write(*,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(4000,name) + goto 9999 + endif + + b_glob => aux_b(:,1) + do i=1, m_problem + b_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_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(4000,name) + goto 9999 + endif + b_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_glob(1:m_problem)) + end if + + call psb_barrier(ictxt) + if (iam==psb_root_) write(*,'("Partition type: block")') + call psb_matdist(aux_A, A, part_block, ictxt, & + & desc_A,b_glob,b,info) + + t2 = psb_wtime() - t1 + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(*,'(" ")') + write(*,'("Time to read and partition matrix : ",es10.4)')t2 + write(*,'(" ")') + end if + +! set RAS with overlap 2 and ILU(0) on the local blocks + + call mld_precinit(P,'AS',info) + call mld_precset(P,mld_sub_ovr_,2,info) + +! build the preconditioner + + t1 = psb_wtime() + + call mld_precbld(A,desc_A,P,info) + + tprec = psb_wtime()-t1 + call psb_amx(ictxt, tprec) + + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_precbld') + goto 9999 + end if + +! set the initial guess + + call psb_geall(x,desc_A,info) + x(:) =0.0 + call psb_geasb(x,desc_A,info) + +! solve Ax=b with preconditioned BiCGSTAB + + call psb_barrier(ictxt) + t1 = psb_wtime() + + call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,istop=2) + + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + call psb_geall(r,desc_A,info) + r(:) =0.0 + call psb_geasb(r,desc_A,info) + call psb_geaxpby(done,b,dzero,r,desc_A,info) + call psb_spmm(-done,A,x,done,r,desc_A,info) + call psb_genrm2s(resmx,r,desc_A,info) + call psb_geamaxs(resmxp,r,desc_A,info) + + amatsize = psb_sizeof(A) + descsize = psb_sizeof(desc_A) + precsize = mld_sizeof(P) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + call psb_sum(ictxt,precsize) + + call mld_precdescr(P,info) + + if (iam==psb_root_) then + write(*,'(" ")') + write(*,'("Matrix: ",A)')mtrx_file + write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Iterations to convergence : ",i6)')iter + write(*,'("Error estimate on exit : ",es10.4)')err + write(*,'("Time to build prec. : ",es10.4)')tprec + write(*,'("Time to solve system : ",es10.4)')t2 + write(*,'("Time per iteration : ",es10.4)')t2/(iter) + write(*,'("Total time : ",es10.4)')t2+tprec + write(*,'("Residual 2-norm : ",es10.4)')resmx + write(*,'("Residual inf-norm : ",es10.4)')resmxp + write(*,'("Total memory occupation for A : ",i10)')amatsize + write(*,'("Total memory occupation for DESC_A : ",i10)')descsize + write(*,'("Total memory occupation for PREC : ",i10)')precsize + end if + + allocate(x_glob(m_problem),r_glob(m_problem),stat=ierr) + if (ierr /= 0) then + write(0,*) 'allocation error: no data collection' + else + call psb_gather(x_glob,x,desc_A,info,root=psb_root_) + call psb_gather(r_glob,r,desc_A,info,root=psb_root_) + if (iam==psb_root_) then + write(0,'(" ")') + write(0,'("Saving x on file")') + write(20,*) 'matrix: ',mtrx_file + write(20,*) 'computed solution on ',np,' processor(s).' + write(20,*) 'iterations to convergence: ',iter + write(20,*) 'error estimate (infinity norm) on exit:', & + & ' ||r||/(||a||||x||+||b||) = ',err + write(20,*) 'max residual = ',resmx, resmxp + write(20,'(a8,4(2x,a20))') 'I','X(I)','R(I)','B(I)' + do i=1,m_problem + write(20,998) i,x_glob(i),r_glob(i),b_glob(i) + enddo + end if + end if +998 format(i8,4(2x,g20.14)) +993 format(i6,4(1x,e12.6)) + +! deallocate the data structures + + call psb_gefree(b, desc_A,info) + call psb_gefree(x, desc_A,info) + call psb_spfree(A, desc_A,info) + call mld_precfree(P,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 + ! + ! get parameters from standard input + ! + subroutine get_parms(ictxt,mtrx,rhs,itmax,tol) + + use psb_base_mod + implicit none + + integer :: ictxt, itmax + real(psb_dpk_) :: tol + character(len=*) :: mtrx, rhs + integer :: iam, np + + call psb_info(ictxt,iam,np) + + if (iam==psb_root_) then + ! read input parameters + call read_data(mtrx,5) + call read_data(rhs,5) + call read_data(itmax,5) + call read_data(tol,5) + end if + + call psb_bcast(ictxt,mtrx) + call psb_bcast(ictxt,rhs) + call psb_bcast(ictxt,itmax) + call psb_bcast(ictxt,tol) + + end subroutine get_parms +end program mld_dexample_ml diff --git a/examples/fileread/mld_dexample_ml.f90 b/examples/fileread/mld_dexample_ml.f90 new file mode 100644 index 00000000..2d7e6327 --- /dev/null +++ b/examples/fileread/mld_dexample_ml.f90 @@ -0,0 +1,348 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.2) +!!$ +!!$ (C) Copyright 2008 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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: mld_dexample_ml.f90 +! +! This sample program solves a linear system by using BiCGStab coupled with +! one of the following multi-level preconditioner, as explained in Section 6.1 +! of the MLD2P4 User's and Reference Guide: +! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) +! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) +! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) +! +! The matrix and the rhs are read from files (if an rhs is not available, the +! unit rhs is set). +! +program mld_dexample_ml + use psb_base_mod + use mld_prec_mod + use psb_krylov_mod + use psb_util_mod + use data_input + + implicit none + +! input file parameters + character(len=40) :: mtrx_file, rhs_file + +! sparse matrices + type(psb_dspmat_type) :: A, aux_A + +! descriptor of sparse matrices + type(psb_desc_type):: desc_A + +! preconditioner + type(mld_dprec_type) :: P + +! right-hand side, solution and residual vectors + real(psb_dpk_), allocatable , save :: b(:), x(:), r(:), & + & x_glob(:), r_glob(:) + real(psb_dpk_), allocatable, target :: aux_b(:,:) + real(psb_dpk_), pointer :: b_glob(:) + +! solver and preconditioner parameters + real(psb_dpk_) :: tol, err + integer :: itmax, iter, istop + integer :: nlev + +! parallel environment parameters + integer :: ictxt, iam, np + +! other variables + integer :: choice + integer :: i,info,j,m_problem,amatsize,descsize,precsize + integer :: ierr, ircode + real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp + character(len=20) :: name + +! initialize the parallel environment + + 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='mld_dexample_ml' + if(psb_get_errstatus() /= 0) goto 9999 + info=0 + call psb_set_errverbosity(2) + +! get parameters + + call get_parms(ictxt,mtrx_file,rhs_file,choice,itmax,tol) + + call psb_barrier(ictxt) + t1 = psb_wtime() + +! read and assemble the matrix A and the right-hand side b +! using PSBLAS routines for sparse matrix / vector management + + if (iam==psb_root_) then + call read_mat(mtrx_file, aux_A, ictxt) + + m_problem = aux_A%m + call psb_bcast(ictxt,m_problem) + + if(rhs_file /= 'NONE') then + ! reading an rhs + call read_rhs(rhs_file,aux_b,ictxt) + end if + + if (psb_size(aux_b,1)==m_problem) then + ! if any rhs were present, broadcast the first one + write(0,'("Ok, got an rhs ")') + b_glob =>aux_b(:,1) + else + write(*,'("Generating an rhs...")') + write(*,'(" ")') + call psb_realloc(m_problem,1,aux_b,ircode) + if (ircode /= 0) then + call psb_errpush(4000,name) + goto 9999 + endif + + b_glob => aux_b(:,1) + do i=1, m_problem + b_glob(i) = 1.d0 + enddo + endif + call psb_bcast(ictxt,b_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(4000,name) + goto 9999 + endif + b_glob =>aux_b(:,1) + call psb_bcast(ictxt,b_glob(1:m_problem)) + end if + + call psb_barrier(ictxt) + if (iam==psb_root_) write(*,'("Partition type: block")') + call psb_matdist(aux_A, A, part_block, ictxt, & + & desc_A,b_glob,b,info) + + t2 = psb_wtime() - t1 + + call psb_amx(ictxt, t2) + + if (iam==psb_root_) then + write(*,'(" ")') + write(*,'("Time to read and partition matrix : ",es10.4)')t2 + write(*,'(" ")') + end if + + select case(choice) + + case(1) + +! initialize the default multi-level preconditioner, i.e. hybrid +! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) +! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU +! on the blocks) as distributed coarse-level solver + + call mld_precinit(P,'ML',info) + + case(2) + +! set a three-level hybrid Schwarz preconditioner, which uses +! block Jacobi (with ILU(0) on the blocks) as post-smoother, +! a coarsest matrix replicated on the processors, and the +! LU factorization from UMFPACK as coarse-level solver + + call mld_precinit(P,'ML',info,nlev=3) + call mld_precset(P,mld_smoother_type_,'BJAC',info) + call mld_precset(P,mld_coarse_mat_,'REPL',info) + call mld_precset(P,mld_coarse_solve_,'UMF',info) + + case(3) + +! set a three-level additive Schwarz preconditioner, which uses +! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and +! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU +! on the blocks) as distributed coarsest-level solver + + call mld_precinit(P,'ML',info,nlev=3) + call mld_precset(P,mld_ml_type_,'ADD',info) + call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info) + call mld_precset(P,mld_coarse_sweeps_,5,info) + + end select + +! build the preconditioner + + call psb_barrier(ictxt) + t1 = psb_wtime() + + call mld_precbld(A,desc_A,P,info) + + tprec = psb_wtime()-t1 + call psb_amx(ictxt, tprec) + + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_precbld') + goto 9999 + end if + +! set the initial guess + + call psb_geall(x,desc_A,info) + x(:) =0.0 + call psb_geasb(x,desc_A,info) + +! solve Ax=b with preconditioned BiCGSTAB + + call psb_barrier(ictxt) + t1 = psb_wtime() + + call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) + + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + call psb_geall(r,desc_A,info) + r(:) =0.0 + call psb_geasb(r,desc_A,info) + call psb_geaxpby(done,b,dzero,r,desc_A,info) + call psb_spmm(-done,A,x,done,r,desc_A,info) + call psb_genrm2s(resmx,r,desc_A,info) + call psb_geamaxs(resmxp,r,desc_A,info) + + amatsize = psb_sizeof(A) + descsize = psb_sizeof(desc_A) + precsize = mld_sizeof(P) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + call psb_sum(ictxt,precsize) + + call mld_precdescr(P,info) + + if (iam==psb_root_) then + write(*,'(" ")') + write(*,'("Matrix: ",A)')mtrx_file + write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Iterations to convergence : ",i6)')iter + write(*,'("Error estimate on exit : ",es10.4)')err + write(*,'("Time to build prec. : ",es10.4)')tprec + write(*,'("Time to solve system : ",es10.4)')t2 + write(*,'("Time per iteration : ",es10.4)')t2/(iter) + write(*,'("Total time : ",es10.4)')t2+tprec + write(*,'("Residual 2-norm : ",es10.4)')resmx + write(*,'("Residual inf-norm : ",es10.4)')resmxp + write(*,'("Total memory occupation for A : ",i10)')amatsize + write(*,'("Total memory occupation for DESC_A : ",i10)')descsize + write(*,'("Total memory occupation for PREC : ",i10)')precsize + end if + + allocate(x_glob(m_problem),r_glob(m_problem),stat=ierr) + if (ierr /= 0) then + write(0,*) 'allocation error: no data collection' + else + call psb_gather(x_glob,x,desc_A,info,root=psb_root_) + call psb_gather(r_glob,r,desc_A,info,root=psb_root_) + if (iam==psb_root_) then + write(0,'(" ")') + write(0,'("Saving x on file")') + write(20,*) 'matrix: ',mtrx_file + write(20,*) 'computed solution on ',np,' processors.' + write(20,*) 'iterations to convergence: ',iter + write(20,*) 'error estimate (infinity norm) on exit:', & + & ' ||r||/(||a||||x||+||b||) = ',err + write(20,*) 'max residual = ',resmx, resmxp + write(20,'(a8,4(2x,a20))') 'I','X(I)','R(I)','B(I)' + do i=1,m_problem + write(20,998) i,x_glob(i),r_glob(i),b_glob(i) + enddo + end if + end if +998 format(i8,4(2x,g20.14)) +993 format(i6,4(1x,e12.6)) + +! deallocate the data structures + + call psb_gefree(b, desc_A,info) + call psb_gefree(x, desc_A,info) + call psb_spfree(A, desc_A,info) + call mld_precfree(P,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 + ! + ! get parameters from standard input + ! + subroutine get_parms(ictxt,mtrx,rhs,choice,itmax,tol) + + use psb_base_mod + implicit none + + integer :: ictxt, choice, itmax + real(psb_dpk_) :: tol + character(len=*) :: mtrx, rhs + integer :: iam, np + + call psb_info(ictxt,iam,np) + + if (iam==psb_root_) then + ! read input parameters + call read_data(mtrx,5) + call read_data(rhs,5) + call read_data(choice,5) + call read_data(itmax,5) + call read_data(tol,5) + end if + + call psb_bcast(ictxt,mtrx) + call psb_bcast(ictxt,rhs) + call psb_bcast(ictxt,choice) + call psb_bcast(ictxt,itmax) + call psb_bcast(ictxt,tol) + + end subroutine get_parms +end program mld_dexample_ml diff --git a/examples/pdegen/Makefile b/examples/pdegen/Makefile new file mode 100644 index 00000000..affefbf8 --- /dev/null +++ b/examples/pdegen/Makefile @@ -0,0 +1,40 @@ +MLDDIR=../.. +include $(MLDDIR)/Make.inc +PSBDIR=$(PSBLASDIR)/lib/ +MLDLIBDIR=$(MLDDIR)/lib +MLD_LIB=-L$(MLDLIBDIR) -lmld_krylov -lmld_prec +PSBLAS_LIB= -L$(PSBDIR) -lpsb_util -lpsb_base +FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBDIR) $(FIFLAG). + +DMOBJS=mld_dexample_ml.o data_input.o +D1OBJS=mld_dexample_1lev.o data_input.o + +EXEDIR=./runs + +all: mld_dexample_ml mld_dexample_1lev + +mld_dexample_ml: $(DMOBJS) + $(F90LINK) $(LINKOPT) $(DMOBJS) -o mld_dexample_ml \ + $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv mld_dexample_ml $(EXEDIR) + +mld_dexample_1lev: $(D1OBJS) + $(F90LINK) $(LINKOPT) $(D1OBJS) -o mld_dexample_1lev \ + $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) + /bin/mv mld_dexample_1lev $(EXEDIR) + +mld_dexample_ml.o: data_input.o +mld_dexample_1lev.o: data_input.o + +.f90.o: + $(MPF90) $(F90COPT) $(FINCLUDES) -c $< + +clean: + /bin/rm -f $(DMOBJS) $(D1OBJS) \ + *$(.mod) $(EXEDIR)/mld_dexample_ml $(EXEDIR)/mld_dexample_1lev + +lib: + (cd ../../; make library) +verycleanlib: + (cd ../../; make veryclean) + diff --git a/examples/pdegen/data_input.f90 b/examples/pdegen/data_input.f90 new file mode 100644 index 00000000..32aae9fb --- /dev/null +++ b/examples/pdegen/data_input.f90 @@ -0,0 +1,91 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.2) +!!$ +!!$ (C) Copyright 2008 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +module data_input + + interface read_data + module procedure read_char, read_int,& + & read_double, read_single + end interface read_data + +contains + + subroutine read_char(val,file) + character(len=*), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),'(a)') val + end subroutine read_char + subroutine read_int(val,file) + integer, intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_int + subroutine read_single(val,file) + use psb_base_mod + real(psb_spk_), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_single + subroutine read_double(val,file) + use psb_base_mod + real(psb_dpk_), intent(out) :: val + integer, intent(in) :: file + character(len=1024) :: charbuf + integer :: idx + read(file,'(a)')charbuf + charbuf = adjustl(charbuf) + idx=index(charbuf,"!") + read(charbuf(1:idx-1),*) val + end subroutine read_double +end module data_input + diff --git a/examples/pdegen/mld_dexample_1lev.f90 b/examples/pdegen/mld_dexample_1lev.f90 new file mode 100644 index 00000000..7f43e870 --- /dev/null +++ b/examples/pdegen/mld_dexample_1lev.f90 @@ -0,0 +1,598 @@ +!!$ +!!$ +!!$ MLD2P4 version 1.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 2.2) +!!$ +!!$ (C) Copyright 2008 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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: mld_dexample_1lev.f90 +! +! This sample program solves a linear system obtained by discretizing a +! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the +! following multi-level preconditioner, as explained in Section 6.1 of +! the MLD2P4 User's and Reference Guide: +! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) +! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) +! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) +! +! The PDE is a general second order equation in 3d +! +! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) +! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! dxdx dydy dzdz dx dy dz +! +! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! +! Example taken from: +! C.T.Kelley +! Iterative Methods for Linear and Nonlinear Equations +! SIAM 1995 +! +! In this sample program the index space of the discretized +! computational domain is first numbered sequentially in a standard way, +! then the corresponding vector is distributed according to a BLOCK +! data distribution. +! +! Boundary conditions are set in a very simple way, by adding +! equations of the form +! +! u(x,y) = exp(-x^2-y^2-z^2) +! +! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. +! +program mld_dexample_1lev + use psb_base_mod + use mld_prec_mod + use psb_krylov_mod + use psb_util_mod + use data_input + + implicit none + + +! sparse matrices + type(psb_dspmat_type) :: A + +! descriptor of sparse matrices + type(psb_desc_type):: desc_A + +! preconditioner + type(mld_dprec_type) :: P + +! right-hand side, solution and residual vectors + real(psb_dpk_), allocatable , save :: b(:), x(:), r(:) + +! solver parameters + real(psb_dpk_) :: tol, err + integer :: itmax, iter, itrace, istop + +! parallel environment parameters + integer :: ictxt, iam, np + +! other variables + integer :: i,info,j,amatsize,descsize,precsize + integer :: idim, nlev, ierr, ircode + real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp + character(len=20) :: name + +! initialize the parallel environment + 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='mld_dexample_ml' + if(psb_get_errstatus() /= 0) goto 9999 + info=0 + call psb_set_errverbosity(2) + +! get parameters + + call get_parms(ictxt,idim,itmax,tol) + +! allocate and fill in the coefficient matrix, rhs and initial guess + + call psb_barrier(ictxt) + t1 = psb_wtime() + call create_matrix(idim,A,b,x,desc_A,part_block,ictxt,info) + t2 = psb_wtime() - t1 + if(info /= 0) then + info=4010 + call psb_errpush(info,name) + goto 9999 + end if + + call psb_amx(ictxt,t2) + if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es10.4)')t2 + if (iam == psb_root_) write(*,'(" ")') + +! set RAS with overlap 2 and ILU(0) on the local blocks + + call mld_precinit(P,'AS',info) + call mld_precset(P,mld_sub_ovr_,2,info) + +! build the preconditioner + + call psb_barrier(ictxt) + t1 = psb_wtime() + + call mld_precbld(A,desc_A,P,info) + + tprec = psb_wtime()-t1 + call psb_amx(ictxt, tprec) + + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_precbld') + goto 9999 + end if + +! set the initial guess + + call psb_geall(x,desc_A,info) + x(:) =0.0 + call psb_geasb(x,desc_A,info) + +! solve Ax=b with preconditioned BiCGSTAB + + call psb_barrier(ictxt) + t1 = psb_wtime() + + call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2) + + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + call psb_geall(r,desc_A,info) + r(:) =0.0 + call psb_geasb(r,desc_A,info) + call psb_geaxpby(done,b,dzero,r,desc_A,info) + call psb_spmm(-done,A,x,done,r,desc_A,info) + call psb_genrm2s(resmx,r,desc_A,info) + call psb_geamaxs(resmxp,r,desc_A,info) + + amatsize = psb_sizeof(A) + descsize = psb_sizeof(desc_A) + precsize = mld_sizeof(P) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + call psb_sum(ictxt,precsize) + + call mld_precdescr(P,info) + + if (iam==psb_root_) then + write(*,'(" ")') + write(*,'("Matrix from PDE example")') + write(*,'("Computed solution on ",i8," processors")')np + write(*,'("Iterations to convergence : ",i6)')iter + write(*,'("Error estimate on exit : ",es10.4)')err + write(*,'("Time to build prec. : ",es10.4)')tprec + write(*,'("Time to solve system : ",es10.4)')t2 + write(*,'("Time per iteration : ",es10.4)')t2/(iter) + write(*,'("Total time : ",es10.4)')t2+tprec + write(*,'("Residual 2-norm : ",es10.4)')resmx + write(*,'("Residual inf-norm : ",es10.4)')resmxp + write(*,'("Total memory occupation for A : ",i10)')amatsize + write(*,'("Total memory occupation for DESC_A : ",i10)')descsize + write(*,'("Total memory occupation for PREC : ",i10)')precsize + end if + + call psb_gefree(b, desc_A,info) + call psb_gefree(x, desc_A,info) + call psb_spfree(A, desc_A,info) + call mld_precfree(P,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 + ! + ! get parameters from standard input + ! + subroutine get_parms(ictxt,idim,itmax,tol) + + use psb_base_mod + implicit none + + integer :: idim, ictxt, itmax + real(psb_dpk_) :: tol + integer :: iam, np + + call psb_info(ictxt,iam,np) + + if (iam==psb_root_) then + ! read input parameters + call read_data(idim,5) + call read_data(itmax,5) + call read_data(tol,5) + end if + + call psb_bcast(ictxt,idim) + call psb_bcast(ictxt,itmax) + call psb_bcast(ictxt,tol) + + end subroutine get_parms + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs + ! + subroutine create_matrix(idim,a,b,xv,desc_a,parts,ictxt,info) + ! + ! Discretize the partial diferential equation + ! + ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) + ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. + ! + ! Boundary conditions are set in a very simple way, by adding + ! equations of the form + ! + ! u(x,y) = exp(-x^2-y^2-z^2) + ! + ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. + ! + use psb_base_mod + implicit none + integer :: idim + integer, parameter :: nbmax=10 + real(psb_dpk_), allocatable :: b(:),xv(:) + type(psb_desc_type) :: desc_a + integer :: ictxt, info + interface + ! .....user passed subroutine..... + subroutine parts(global_indx,n,np,pv,nv) + implicit none + integer, intent(in) :: global_indx, n, np + integer, intent(out) :: nv + integer, intent(out) :: pv(*) + end subroutine parts + end interface + ! local variables + type(psb_dspmat_type) :: a + real(psb_dpk_) :: zt(nbmax),glob_x,glob_y,glob_z + integer :: m,n,nnz,glob_row + integer :: x,y,z,ia,indx_owner + integer :: np, iam + integer :: element + integer :: nv, inv + integer, allocatable :: irow(:),icol(:) + real(psb_dpk_), allocatable :: val(:) + integer, allocatable :: prv(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_dpk_) :: deltah + real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t1, t2, t3, tins, tasb + real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 + external :: a1, a2, a3, a4, b1, b2, b3 + integer :: err_act + + character(len=20) :: name + + info = 0 + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + deltah = 1.d0/(idim-1) + + ! initialize array descriptor and sparse matrix storage; provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(0,'("Generating Matrix (size=",i0x,")...")')n + + call psb_cdall(ictxt,desc_a,info,mg=n,parts=parts) + call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + call psb_geall(b,desc_a,info) + call psb_geall(xv,desc_a,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nbmax),irow(20*nbmax),& + &icol(20*nbmax),prv(np),stat=info) + if (info /= 0 ) then + info=4000 + call psb_errpush(info,name) + goto 9999 + endif + + tins = 0.d0 + call psb_barrier(ictxt) + t1 = psb_wtime() + + ! loop over rows belonging to current process in a block + ! distribution. + + ! icol(1)=1 + do glob_row = 1, n + call parts(glob_row,n,np,prv,nv) + do inv = 1, nv + indx_owner = prv(inv) + if (indx_owner == iam) then + ! local matrix pointer + element=1 + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + x = glob_row/(idim*idim) + else + x = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then + y = (glob_row-(x-1)*idim*idim)/idim + else + y = (glob_row-(x-1)*idim*idim)/idim+1 + endif + z = glob_row-(x-1)*idim*idim-(y-1)*idim + ! glob_x, glob_y, glob_x coordinates + glob_x=x*deltah + glob_y=y*deltah + glob_z=z*deltah + + ! check on boundary points + zt(1) = 0.d0 + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + if (x==1) then + val(element)=-b1(glob_x,glob_y,glob_z)& + & -a1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element)) + else + val(element)=-b1(glob_x,glob_y,glob_z)& + & -a1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-2)*idim*idim+(y-1)*idim+(z) + element=element+1 + endif + ! term depending on (x,y-1,z) + if (y==1) then + val(element)=-b2(glob_x,glob_y,glob_z)& + & -a2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b2(glob_x,glob_y,glob_z)& + & -a2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-1)*idim*idim+(y-2)*idim+(z) + element=element+1 + endif + ! term depending on (x,y,z-1) + if (z==1) then + val(element)=-b3(glob_x,glob_y,glob_z)& + & -a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b3(glob_x,glob_y,glob_z)& + & -a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1) + element=element+1 + endif + ! term depending on (x,y,z) + val(element)=2*b1(glob_x,glob_y,glob_z)& + & +2*b2(glob_x,glob_y,glob_z)& + & +2*b3(glob_x,glob_y,glob_z)& + & +a1(glob_x,glob_y,glob_z)& + & +a2(glob_x,glob_y,glob_z)& + & +a3(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-1)*idim*idim+(y-1)*idim+(z) + element=element+1 + ! term depending on (x,y,z+1) + if (z==idim) then + val(element)=-b1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b1(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1) + element=element+1 + endif + ! term depending on (x,y+1,z) + if (y==idim) then + val(element)=-b2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) + else + val(element)=-b2(glob_x,glob_y,glob_z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element)=(x-1)*idim*idim+(y)*idim+(z) + element=element+1 + endif + ! term depending on (x+1,y,z) + if (x=0.0) then - call mld_precset(pre,mld_aggr_damp_,precs(pp)%omega,info) - end if - call mld_precset(pre,mld_ml_type_, precs(pp)%mltype, info) - call mld_precset(pre,mld_aggr_alg_, precs(pp)%aggr, info) - call mld_precset(pre,mld_aggr_kind_, precs(pp)%smthkind, info) - call mld_precset(pre,mld_smoother_pos_, precs(pp)%smthpos, info) - call mld_precset(pre,mld_coarse_solve_, precs(pp)%cslv, info) - call mld_precset(pre,mld_coarse_subsolve_, precs(pp)%ftype2, info) - call mld_precset(pre,mld_coarse_fillin_, precs(pp)%fill2, info) - call mld_precset(pre,mld_coarse_fthrs_, precs(pp)%thr2, info) - call mld_precset(pre,mld_coarse_sweeps_, precs(pp)%jswp, info) - call mld_precset(pre,mld_coarse_mat_, precs(pp)%cmat, info) - else - call mld_precinit(pre,precs(pp)%lv1,info) - end if - call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype1, info) - call mld_precset(pre,mld_sub_restr_, precs(pp)%restr, info) - call mld_precset(pre,mld_sub_prol_, precs(pp)%prol, info) - call mld_precset(pre,mld_sub_ovr_, precs(pp)%novr, info) - call mld_precset(pre,mld_sub_fillin_, precs(pp)%fill1, info) - call mld_precset(pre,mld_fact_thrs_, precs(pp)%thr1, info) - - - ! setting initial guess to zero - call psb_geaxpby(dzero,b_col,dzero,x_col,desc_a,info) - - ! building the preconditioner -!!$ write(0,*) 'Check in main program on hasv in:',allocated(desc_a%hashv) -!!$ call flush(0) - call psb_barrier(ictxt) - t1 = psb_wtime() - call mld_precbld(a,desc_a,pre,info) - tprec = psb_wtime()-t1 - if (info /= 0) then - write(0,*) 'INFO from precbld ',info - call psb_error() - goto 9999 - end if - if (psb_get_errstatus() /= 0) then - write(0,*) 'INFO from precbld ',info - call psb_error() - goto 9999 - end if - -!!$ write(0,*) 'Check in main program on hasv out:',allocated(desc_a%hashv) -!!$ call flush(0) - call psb_amx(ictxt,tprec) -!!$ call psb_csprt(40+iam,a,head='% (A)') -!!$ if (iam == psb_root_) then -!!$ write(*,'("Matrix : ",a)') mtrx(nm) -!!$ write(*,'("RHS : ",a)') rhs(nm) -!!$ write(*,'("Method : ",a)') kmethd -!!$ write(*,'("Preconditioner : ",a)') precs(pp)%descr -!!$ call mld_precdescr(pre) -!!$ call flush(6) -!!$ end if - iparm = 0 - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_krylov(kmethd,a,pre,b_col,x_col,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,& - & irst=irst,istop=istopc) - call psb_barrier(ictxt) - t2 = psb_wtime() - t1 - call psb_amx(ictxt,t2) - if (info/=0) then - write(0,*) 'INFO from solver ',info - call psb_errpush(4010,name,a_err='psb_SOLVER') - goto 9999 - end if -!!$ write(0,*) iam,'Done Solver' -!!$ call flush(0) -!!$ call flush(6) -!!$ call psb_barrier(ictxt) - if(iam == psb_root_.and.out2) & - & write(10,'(a20,2(1x,i3),1x,i5,3(1x,g9.4),1x,a8,1x,a)') & - & mtrx(nm),np,precs(pp)%novr,iter,tprec,t2,t2+tprec,& - & trim(kmethd),trim(precs(pp)%descr) - if(iam == psb_root_) & - & write(0,'(a20,2(1x,i3),1x,i5,3(1x,g9.4),1x,a8,1x,a)') & - & mtrx(nm),np,precs(pp)%novr,iter,tprec,t2,t2+tprec,& - & trim(kmethd),trim(precs(pp)%descr) - if (nt.lt.ntry) call mld_precfree(pre,info) - if((t2+tprec).lt.mttot) then - mtslv=t2 - mtprec=tprec - mttot=t2+tprec - end if - end do -!!$ call psb_csprt(50+iam,a,head='% (A)') -!!$ write(0,*) 'Check hashv after precfree:',allocated(desc_a%hashv) -!!$ call flush(0) - - call psb_geaxpby(done,b_col,dzero,r_col,desc_a,info) - call psb_spmm(-done,a,x_col,done,r_col,desc_a,info) - call psb_genrm2s(resmx,r_col,desc_a,info) - call psb_geamaxs(resmxp,r_col,desc_a,info) - - ipsize = mld_sizeof(pre) - asize = psb_sizeof(a) - cdsize = psb_sizeof(desc_a) - call psb_sum(ictxt,ipsize) - call psb_sum(ictxt,asize) - call psb_sum(ictxt,cdsize) - - if (iam == psb_root_) then - write(*,'("Matrix : ",a)') mtrx(nm) - write(*,'("RHS : ",a)') rhs(nm) - write(*,'("Method : ",a)') kmethd - write(*,'("Preconditioner : ",a)') precs(pp)%descr - call mld_precdescr(pre) - write(*,'("Computed solution on ",i4," processors")')np - write(*,'(" ")') - write(*,'("Iterations to convergence: ",i6)') iter - write(*,'("Error indicator on exit : ",g9.4)') err - write(*,'("Time to buil prec. : ",es10.4)')mtprec - write(*,'("Time to solve matrix : ",es10.4)')mtslv - write(*,'("Time per iteration : ",es10.4)')mtslv/(iter) - write(*,'("Total time : ",es10.4)')mttot - write(*,'("Residual norm 2 : ",es10.4)')resmx - write(*,'("Residual norm inf : ",es10.4)')resmxp - write(*,'("Total memory occupation for A: ",i10)')asize - write(*,'("Total memory occupation for DESC_A: ",i10)')cdsize - write(*,'("Total memory occupation for PRE: ",i10)')ipsize - - write(*,'(72("="))') - write(*,'(" ")') - write(*,'(" ")') - write(*,'(" ")') - - if(out1) write(8,'(a20,2(1x,i3),1x,i5,5(1x,g9.4),1x,a8,1x,a)') mtrx(nm),& - & np,precs(pp)%novr,& - & iter,mtprec,mtslv,mttot,resmx,resmxp,& - & trim(kmethd),trim(precs(pp)%descr) - end if - - call mld_precfree(pre,info) - - if (.false.) then - allocate(x_col_glob(m_problem),r_col_glob(m_problem),stat=ierr) - if (ierr /= 0) then - write(0,*) 'allocation error: no data collection' - else - call psb_gather(x_col_glob,x_col,desc_a,info,root=psb_root_) - call psb_gather(r_col_glob,r_col,desc_a,info,root=psb_root_) - if (iam == psb_root_) then - write(0,'(" ")') - write(0,'("Saving x on file")') - write(outf3,'(a,a,a)')trim(mtrx(nm)),'.psb_sol.',& - & psb_tolower(trim(precs(pp)%descr)) - open(20,file=outf3) - write(20,*) 'matrix: ',mtrx(nm) - write(20,*) 'computed solution on ',np,' processors.' - write(20,*) 'iterations to convergence: ',iter - write(20,*) 'error indicator (infinity norm) on exit:', & - & ' ||r||/(||a||||x||+||b||) = ',err - write(20,*) 'max residual = ',resmx, resmxp - do i=1,m_problem - write(20,998) i,x_col_glob(i),b_col_glob(i) - enddo - close(20) - end if - end if - end if -998 format(i8,4(2x,g20.14)) -993 format(i6,4(1x,e12.6)) - - -!!$ if (pp == 1) call psb_csprt(40+iam,a,head='% (A)') - end do - - 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) - deallocate(r_col,stat=info) - deallocate(aux_b,stat=info) - if (ipart==0) then - deallocate(ivg,ipv,stat=info) - endif - end do - -9999 continue - if(info /= 0) then - call psb_error(ictxt) - end if - call psb_exit(ictxt) - stop - - -contains - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(icontxt,irst,irnum,ntry,nmat,mtrx,rhs,kmethd,nprecs,precs,ipart,& - & afmt,istopc,itmax,itrace,eps,outf1,outf2) - - use psb_base_mod - implicit none - - integer :: icontxt - character(len=20) :: kmethd - character(len=80) :: outf1, outf2 - character(len=20),allocatable :: mtrx(:), rhs(:) - type(precdata),allocatable :: precs(:) - integer :: iret, istopc,itmax,itrace,ipart,nmat,nprecs,irst,irnum,ntry - character(len=1024) :: charbuf - real(psb_dpk_) :: eps, omega,thr1,thr2 - character :: afmt*5, lv1*10, lv2*10, pdescr*40 - integer :: iam, nm, np, i, idx - integer, parameter :: npparms=15 - integer :: inparms(40), ip, pparms(npparms) - - call psb_info(icontxt,iam,np) - - if (iam==psb_root_) then - ! read input parameters - read(*,*) outf1 - read(*,*) outf2 - read(*,*) kmethd - read(*,*) eps - read(*,*) afmt - - call psb_bcast(icontxt,kmethd) - call psb_bcast(icontxt,eps,0) - - call psb_bcast(icontxt,afmt) - - read(*,*) ipart - read(*,*) itmax - read(*,*) itrace - read(*,*) istopc - read(*,*) irst - read(*,*) irnum - read(*,*) ntry - read(*,*) nprecs - ! broadcast parameters to all processors - - inparms(1) = ipart - inparms(2) = itmax - inparms(3) = itrace - inparms(4) = istopc - inparms(5) = irst - inparms(6) = irnum - inparms(7) = ntry - call psb_bcast(icontxt,inparms(1:7)) - - call psb_bcast(icontxt,nprecs) - - allocate(precs(nprecs)) - - do np=1,nprecs - read(*,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf," ") - read(charbuf(1:idx-1),'(a)')lv1 - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx-1),'(a)')lv2 - charbuf=adjustl(charbuf(idx:)) - do i=1, npparms - idx=index(charbuf," ") - read(charbuf(1:idx),*) pparms(i) - charbuf=adjustl(charbuf(idx:)) - end do - idx=index(charbuf," ") - read(charbuf(1:idx),*) omega - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx),*) thr1 - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx),*) thr2 - charbuf=adjustl(charbuf(idx:)) - read(charbuf,'(a)') pdescr - - call psb_bcast(icontxt,pdescr) - precs(np)%descr=pdescr - - call psb_bcast(icontxt,lv1) - call psb_bcast(icontxt,lv2) - call psb_bcast(icontxt,pparms(1:npparms)) - call psb_bcast(icontxt,omega) - call psb_bcast(icontxt,thr1) - call psb_bcast(icontxt,thr2) - - precs(np)%lv1 = lv1 - precs(np)%lv2 = lv2 - precs(np)%novr = pparms(1) - precs(np)%restr = pparms(2) - precs(np)%prol = pparms(3) - precs(np)%ftype1 = pparms(4) - precs(np)%fill1 = pparms(5) - precs(np)%mltype = pparms(6) - precs(np)%aggr = pparms(7) - precs(np)%smthkind = pparms(8) - precs(np)%cmat = pparms(9) - precs(np)%smthpos = pparms(10) - precs(np)%cslv = pparms(11) - precs(np)%ftype2 = pparms(12) - precs(np)%fill2 = pparms(13) - precs(np)%jswp = pparms(14) - precs(np)%nlev = pparms(15) - precs(np)%omega = omega - precs(np)%thr1 = thr1 - precs(np)%thr2 = thr2 - end do - - read(*,*) nmat - call psb_bcast(icontxt,nmat) - allocate(mtrx(nmat),rhs(nmat)) - - do nm=1, nmat - read(*,'(a)') charbuf - charbuf=adjustl(charbuf) - idx=index(charbuf," ") - mtrx(nm)=charbuf(1:idx-1) - rhs(nm)=adjustl(charbuf(idx+1:)) - call psb_bcast(icontxt,mtrx(nm)) - call psb_bcast(icontxt,rhs(nm)) - end do - - else - ! receive parameters - call psb_bcast(icontxt,kmethd) - call psb_bcast(icontxt,eps) - - call psb_bcast(icontxt,afmt) - - call psb_bcast(icontxt,inparms(1:7)) - - ipart = inparms(1) - itmax = inparms(2) - itrace = inparms(3) - istopc = inparms(4) - irst = inparms(5) - irnum = inparms(6) - ntry = inparms(7) - - call psb_bcast(icontxt,nprecs) - allocate(precs(nprecs)) - - do np=1,nprecs - call psb_bcast(icontxt,pdescr) - precs(np)%descr=pdescr - - call psb_bcast(icontxt,lv1) - - call psb_bcast(icontxt,lv2) - - call psb_bcast(icontxt,pparms(1:npparms)) - call psb_bcast(icontxt,omega) - call psb_bcast(icontxt,thr1) - call psb_bcast(icontxt,thr2) - - precs(np)%lv1 = lv1 - precs(np)%lv2 = lv2 - precs(np)%novr = pparms(1) - precs(np)%restr = pparms(2) - precs(np)%prol = pparms(3) - precs(np)%ftype1 = pparms(4) - precs(np)%fill1 = pparms(5) - precs(np)%mltype = pparms(6) - precs(np)%aggr = pparms(7) - precs(np)%smthkind = pparms(8) - precs(np)%cmat = pparms(9) - precs(np)%smthpos = pparms(10) - precs(np)%cslv = pparms(11) - precs(np)%ftype2 = pparms(12) - precs(np)%fill2 = pparms(13) - precs(np)%jswp = pparms(14) - precs(np)%nlev = pparms(15) - precs(np)%omega = omega - precs(np)%thr1 = thr1 - precs(np)%thr2 = thr2 - end do - - - call psb_bcast(icontxt,nmat) - allocate(mtrx(nmat),rhs(nmat)) - - do nm=1,nmat - - call psb_bcast(icontxt,mtrx(nm)) - call psb_bcast(icontxt,rhs(nm)) - - end do - - end if - - end subroutine get_parms - subroutine pr_usage(iout) - integer iout - write(iout, *) ' number of parameters is incorrect!' - write(iout, *) ' use: hb_sample mtrx_file methd prec [ptype & - &itmax istopc itrace]' - write(iout, *) ' where:' - write(iout, *) ' mtrx_file is stored in hb format' - write(iout, *) ' methd may be: cgstab ' - write(iout, *) ' itmax max iterations [500] ' - write(iout, *) ' istopc stopping criterion [1] ' - write(iout, *) ' itrace 0 (no tracing, default) or ' - write(iout, *) ' >= 0 do tracing every itrace' - write(iout, *) ' iterations ' - write(iout, *) ' prec may be: ilu diagsc none' - write(iout, *) ' ptype partition strategy default 0' - write(iout, *) ' 0: block partition ' - end subroutine pr_usage - -end program df_bench - - - - - diff --git a/test/fileread/runs/drt.sh b/test/fileread/runs/drt.sh deleted file mode 100755 index caeab16c..00000000 --- a/test/fileread/runs/drt.sh +++ /dev/null @@ -1,84 +0,0 @@ -#!/bin/sh - -# -# -# NUmber of attempts for each configuration -ntry=1 - -export GFORTRAN_UNBUFFERED_ALL=y -date=`date +%Y%m%d%H%M%S` - - - -for np in $* -do - -# 3rd batch: 14bis,64bis 4 sweeps -echo "mpirun -np $np -machinefile locm df_bench >>log.part$part.ren$renum.${np}p" -#/usr/local/mpich-gcc42/bin/mpirun -np $np -machinefile locm df_bench <>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <>dat.out.kiva.$date -cat stat.${np}p.$date >>dat.stat.kiva.$date - -done - - diff --git a/test/fileread/zf_bench.f90 b/test/fileread/zf_bench.f90 deleted file mode 100644 index 17de70ff..00000000 --- a/test/fileread/zf_bench.f90 +++ /dev/null @@ -1,594 +0,0 @@ -program zf_bench - use psb_base_mod - use psb_util_mod - use mld_prec_mod - use psb_krylov_mod - implicit none - - ! input parameters - character(len=20) :: kmethd - character(len=80) :: outf1, outf2, outf3 - character(len=20), allocatable :: mtrx(:),rhs(:) - type precdata - character(len=10) :: lv1, lv2 ! First and second level prec type - integer :: nlev ! - integer :: novr ! number of overlapping levels - integer :: restr ! restriction over application of as - integer :: prol ! prolongation over application of as - integer :: ftype1 ! Factorization type: ILU, SuperLU, UMFPACK. - integer :: fill1 ! Fill-in for factorization 1 - real(psb_dpk_) :: thr1 ! Threshold for fact. 1 ILU(T) - integer :: mltype ! additive or multiplicative 2nd level prec - integer :: aggr ! local or global aggregation - integer :: smthkind ! smoothing type - integer :: cmat ! coarse mat - integer :: smthpos ! pre, post, both smoothing - integer :: glbsmth ! global smoothing - integer :: cslv ! Coarse solver: BJAC, SuperLU_Dist, UMF. - integer :: ftype2 ! Factorization type: ILU, SuperLU, UMFPACK. - integer :: fill2 ! Fill-in for factorization 1 - real(psb_dpk_) :: thr2 ! Threshold for fact. 1 ILU(T) - integer :: jswp ! Jacobi sweeps - real(psb_dpk_) :: omega ! smoother omega - real(psb_dpk_) :: athres ! smoother aggregation - character(len=40) :: descr ! verbose description of the prec - end type precdata - type(precdata), allocatable :: precs(:) - - - ! sparse matrices - type(psb_zspmat_type) :: a, aux_a - - ! preconditioner data - type(mld_zprec_type) :: pre - integer :: igsmth, matop, novr - - ! dense matrices - complex(psb_dpk_), allocatable, target :: aux_b(:,:), d(:) - complex(psb_dpk_), allocatable , save :: b_col(:), x_col(:), r_col(:), & - & x_col_glob(:), r_col_glob(:) - complex(psb_dpk_), pointer :: b_col_glob(:) - - ! communications data structure - type(psb_desc_type):: desc_a - - ! blacs variables - integer :: ictxt, iam, np - logical :: out1, out2 - - ! solver paramters - integer :: iter, itmax, ierr, itrace, ircode, ipart,nlev,& - & methd, istopc, iprec, ml, irnum, irst, ntry, nmat, ilev,ipsize,asize,cdsize - real(psb_dpk_) :: err, eps - - character(len=5) :: afmt - character(len=20) :: name - integer :: iparm(20) - - ! other variables - integer :: i,info,j,m_problem, nm, nt - integer :: internal, m,ii,nnzero, nprecs, pp - real(psb_dpk_) :: t1, t2, tprec, r_amax, b_amax,& - &scale,resmx,resmxp, mttot, mtslv, mtprec - integer :: nrhs, nrow, n_row, dim, nv, ne - integer, allocatable :: ipv(:), neigh(:), ivg(:) - - - 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='zf_bench' - if(psb_get_errstatus() /= 0) goto 9999 - info=0 - call psb_set_errverbosity(2) - ! - ! get parameters - ! - call get_parms(ictxt,irst,irnum,ntry,nmat,mtrx,rhs,kmethd,nprecs,precs,& - & ipart,afmt,istopc,itmax,itrace,eps,outf1,outf2) - - if(iam == psb_root_) then - if(outf1 /= 'NONE') then - open(8,file=outf1,action='write') - out1=.true. - else - out1=.false. - end if - - if(outf2 /= 'NONE') then - open(10,file=outf2,action='write') - out2=.true. - else - out2=.false. - end if - end if - - do nm=1, nmat - - if(iam == psb_root_) write(*,'(25("=")," ",a20," ",25("="))')mtrx(nm) - 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 - call read_mat(mtrx(nm), aux_a, ictxt) - - m_problem = aux_a%m - call psb_bcast(ictxt,m_problem) - - if(rhs(nm) /= 'none') then - ! reading an rhs - call read_rhs(rhs(nm),aux_b,ictxt) - end if - - if (allocated(aux_b).and.size(aux_b,1)==m_problem) then - ! if any rhs were present, broadcast the first one - b_col_glob =>aux_b(:,1) - else - allocate(aux_b(m_problem,1), stat=ircode) - if (ircode /= 0) then - call psb_errpush(4000,name) - goto 9999 - endif - - b_col_glob => aux_b(:,1) - 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) - allocate(aux_b(m_problem,1), stat=ircode) - if (ircode /= 0) then - call psb_errpush(4000,name) - goto 9999 - endif - b_col_glob =>aux_b(:,1) - call 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(*,'("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, ivg, ictxt, & - & desc_a,b_col_glob,b_col,info,fmt=afmt) - else if (ipart == 2) then - if (iam == psb_root_) then - call build_mtpart(aux_a%m,aux_a%fida,aux_a%ia1,aux_a%ia2,np) - endif - call psb_barrier(ictxt) - call distr_mtpart(psb_root_,ictxt) - call getv_mtpart(ivg) - call psb_matdist(aux_a, a, ivg, ictxt, & - & desc_a,b_col_glob,b_col,info,fmt=afmt) - call free_part(info) - else - call psb_matdist(aux_a, a, part_block, ictxt, & - & desc_a,b_col_glob,b_col,info,fmt=afmt) - end if - - call psb_geall(x_col,desc_a,info) - x_col(:) =0.0 - call psb_geasb(x_col,desc_a,info) - call psb_geall(r_col,desc_a,info) - r_col(:) =0.0 - call psb_geasb(r_col,desc_a,info) - t2 = psb_wtime() - t1 - -!!$ call psb_csprt(10+iam,a,head='% (A)') - call psb_sp_free(aux_a,info) - - call psb_amx(ictxt, t2) - -!!$ call psb_csprt(20+iam,a,head='% (A)') - ! - ! prepare the preconditioning matrix. note the availability - ! of optional parameters - ! - - - do pp=1, nprecs - - mttot=1.d300 - - do nt=1,ntry - - if (precs(pp)%lv2(1:2) == 'ml') then - if (precs(pp)%nlev < 2) then - write(0,*) 'Inconsistent number of levels ',precs(pp)%nlev,& - & ' forcing 2' - precs(pp)%nlev = 2 - end if - nlev = precs(pp)%nlev - call mld_precinit(pre,precs(pp)%lv2,info,nlev=nlev) - ! Defaults are OK for all intermediate levels. Only fix last level. - if (precs(pp)%omega>=0.0) then - call mld_precset(pre,mld_aggr_damp_,precs(pp)%omega,info) - end if - call mld_precset(pre,mld_ml_type_, precs(pp)%mltype, info) - call mld_precset(pre,mld_aggr_alg_, precs(pp)%aggr, info) - call mld_precset(pre,mld_aggr_kind_, precs(pp)%smthkind, info) - call mld_precset(pre,mld_smoother_pos_, precs(pp)%smthpos, info) - call mld_precset(pre,mld_coarse_solve_, precs(pp)%cslv, info) - call mld_precset(pre,mld_coarse_subsolve_, precs(pp)%ftype2, info) - call mld_precset(pre,mld_coarse_fillin_, precs(pp)%fill2, info) - call mld_precset(pre,mld_coarse_fthrs_, precs(pp)%thr2, info) - call mld_precset(pre,mld_coarse_sweeps_, precs(pp)%jswp, info) - call mld_precset(pre,mld_coarse_mat_, precs(pp)%cmat, info) - else - call mld_precinit(pre,precs(pp)%lv1,info) - end if - call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype1, info) - call mld_precset(pre,mld_sub_restr_, precs(pp)%restr, info) - call mld_precset(pre,mld_sub_prol_, precs(pp)%prol, info) - call mld_precset(pre,mld_sub_ovr_, precs(pp)%novr, info) - call mld_precset(pre,mld_sub_fillin_, precs(pp)%fill1, info) - call mld_precset(pre,mld_fact_thrs_, precs(pp)%thr1, info) - - - ! setting initial guess to zero - call psb_geaxpby(zzero,b_col,zzero,x_col,desc_a,info) - - ! building the preconditioner - call psb_barrier(ictxt) - t1 = psb_wtime() - call mld_precbld(a,desc_a,pre,info) - tprec = psb_wtime()-t1 - if (info /= 0) then - write(0,*) 'INFO from precbld ',info - call psb_error() - goto 9999 - end if - - - call psb_amx(ictxt,tprec) -!!$ call psb_csprt(40+iam,a,head='% (A)') -!!$ if (iam == psb_root_) then -!!$ write(*,'("Matrix : ",a)') mtrx(nm) -!!$ write(*,'("RHS : ",a)') rhs(nm) -!!$ write(*,'("Method : ",a)') kmethd -!!$ write(*,'("Preconditioner : ",a)') precs(pp)%descr -!!$ call mld_precdescr(pre) -!!$ end if - iparm = 0 - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_krylov(kmethd,a,pre,b_col,x_col,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,& - & irst=irst,istop=istopc) - call psb_barrier(ictxt) - t2 = psb_wtime() - t1 - call psb_amx(ictxt,t2) - if (info/=0) then - write(0,*) 'INFO from solver ',info - call psb_errpush(4010,name,a_err='psb_SOLVER') - goto 9999 - end if -!!$ write(0,*) iam,'Done Solver' -!!$ call flush(0) -!!$ call flush(6) -!!$ call psb_barrier(ictxt) - if(iam == psb_root_.and.out2) & - & write(10,'(a20,2(1x,i3),1x,i5,3(1x,g9.4),1x,a8,1x,a)') & - & mtrx(nm),np,precs(pp)%novr,iter,tprec,t2,t2+tprec,& - & trim(kmethd),trim(precs(pp)%descr) - if(iam == psb_root_) & - & write(0,'(a20,2(1x,i3),1x,i5,3(1x,g9.4),1x,a8,1x,a)') & - & mtrx(nm),np,precs(pp)%novr,iter,tprec,t2,t2+tprec,& - & trim(kmethd),trim(precs(pp)%descr) - if (nt < ntry) call mld_precfree(pre,info) - if((t2+tprec) < mttot) then - mtslv=t2 - mtprec=tprec - mttot=t2+tprec - end if - - end do -!!$ call psb_csprt(50+iam,a,head='% (A)') - - call psb_geaxpby(zone,b_col,zzero,r_col,desc_a,info) - call psb_spmm(-zone,a,x_col,zone,r_col,desc_a,info) - call psb_genrm2s(resmx,r_col,desc_a,info) - call psb_geamaxs(resmxp,r_col,desc_a,info) - - ipsize = mld_sizeof(pre) - asize = psb_sizeof(a) - cdsize = psb_sizeof(desc_a) - call psb_sum(ictxt,ipsize) - call psb_sum(ictxt,asize) - call psb_sum(ictxt,cdsize) - - if (iam == psb_root_) then - write(*,'("Matrix : ",a)') mtrx(nm) - write(*,'("RHS : ",a)') rhs(nm) - write(*,'("Method : ",a)') kmethd - write(*,'("Preconditioner : ",a)') precs(pp)%descr - call mld_precdescr(pre) - write(*,'("Computed solution on ",i4," processors")')np - write(*,'(" ")') - write(*,'("Iterations to convergence: ",i6)') iter - write(*,'("Error indicator on exit : ",g9.4)') err - write(*,'("Time to buil prec. : ",es10.4)')mtprec - write(*,'("Time to solve matrix : ",es10.4)')mtslv - write(*,'("Time per iteration : ",es10.4)')mtslv/(iter) - write(*,'("Total time : ",es10.4)')mttot - write(*,'("Residual norm 2 : ",es10.4)')resmx - write(*,'("Residual norm inf : ",es10.4)')resmxp - write(*,'("Total memory occupation for A: ",i10)')asize - write(*,'("Total memory occupation for DESC_A: ",i10)')cdsize - write(*,'("Total memory occupation for PRE: ",i10)')ipsize - - write(*,'(72("="))') - write(*,'(" ")') - write(*,'(" ")') - write(*,'(" ")') - - if(out1) write(8,'(a20,2(1x,i3),1x,i5,5(1x,g9.4),1x,a8,1x,a)') mtrx(nm),& - & np,precs(pp)%novr,& - & iter,mtprec,mtslv,mttot,resmx,resmxp,& - & trim(kmethd),trim(precs(pp)%descr) - end if - - call mld_precfree(pre,info) - - end do - - 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) - deallocate(r_col,stat=info) - deallocate(aux_b,stat=info) - if (ipart==0) then - deallocate(ivg,ipv,stat=info) - endif - end do - -9999 continue - if(info /= 0) then - call psb_error(ictxt) - end if - call psb_exit(ictxt) - stop - -contains - - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(icontxt,irst,irnum,ntry,nmat,mtrx,rhs,kmethd,nprecs,precs,ipart,& - & afmt,istopc,itmax,itrace,eps,outf1,outf2) - - use psb_base_mod - implicit none - - integer :: icontxt - character(len=20) :: kmethd - character(len=80) :: outf1, outf2 - character(len=20),allocatable :: mtrx(:), rhs(:) - type(precdata),allocatable :: precs(:) - integer :: iret, istopc,itmax,itrace,ipart,nmat,nprecs,irst,irnum,ntry - character(len=1024) :: charbuf - real(psb_dpk_) :: eps, omega,thr1,thr2 - character :: afmt*5, lv1*10, lv2*10, pdescr*40 - integer :: iam, nm, np, i, idx - integer, parameter :: npparms=15 - integer :: inparms(40), ip, pparms(npparms) - - call psb_info(icontxt,iam,np) - - if (iam==psb_root_) then - ! read input parameters - read(*,*) outf1 - read(*,*) outf2 - read(*,*) kmethd - read(*,*) eps - read(*,*) afmt - - call psb_bcast(icontxt,kmethd) - call psb_bcast(icontxt,eps,0) - - call psb_bcast(icontxt,afmt) - - read(*,*) ipart - read(*,*) itmax - read(*,*) itrace - read(*,*) istopc - read(*,*) irst - read(*,*) irnum - read(*,*) ntry - read(*,*) nprecs - ! broadcast parameters to all processors - - inparms(1) = ipart - inparms(2) = itmax - inparms(3) = itrace - inparms(4) = istopc - inparms(5) = irst - inparms(6) = irnum - inparms(7) = ntry - call psb_bcast(icontxt,inparms(1:7)) - - call psb_bcast(icontxt,nprecs) - - allocate(precs(nprecs)) - - do np=1,nprecs - read(*,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf," ") - read(charbuf(1:idx-1),'(a)')lv1 - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx-1),'(a)')lv2 - charbuf=adjustl(charbuf(idx:)) - do i=1, npparms - idx=index(charbuf," ") - read(charbuf(1:idx),*) pparms(i) - charbuf=adjustl(charbuf(idx:)) - end do - idx=index(charbuf," ") - read(charbuf(1:idx),*) omega - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx),*) thr1 - charbuf=adjustl(charbuf(idx:)) - idx=index(charbuf," ") - read(charbuf(1:idx),*) thr2 - charbuf=adjustl(charbuf(idx:)) - read(charbuf,'(a)') pdescr - - call psb_bcast(icontxt,pdescr) - precs(np)%descr=pdescr - - call psb_bcast(icontxt,lv1) - call psb_bcast(icontxt,lv2) - call psb_bcast(icontxt,pparms(1:npparms)) - call psb_bcast(icontxt,omega) - call psb_bcast(icontxt,thr1) - call psb_bcast(icontxt,thr2) - - precs(np)%lv1 = lv1 - precs(np)%lv2 = lv2 - precs(np)%novr = pparms(1) - precs(np)%restr = pparms(2) - precs(np)%prol = pparms(3) - precs(np)%ftype1 = pparms(4) - precs(np)%fill1 = pparms(5) - precs(np)%mltype = pparms(6) - precs(np)%aggr = pparms(7) - precs(np)%smthkind = pparms(8) - precs(np)%cmat = pparms(9) - precs(np)%smthpos = pparms(10) - precs(np)%cslv = pparms(11) - precs(np)%ftype2 = pparms(12) - precs(np)%fill2 = pparms(13) - precs(np)%jswp = pparms(14) - precs(np)%nlev = pparms(15) - precs(np)%omega = omega - precs(np)%thr1 = thr1 - precs(np)%thr2 = thr2 - end do - - read(*,*) nmat - call psb_bcast(icontxt,nmat) - allocate(mtrx(nmat),rhs(nmat)) - - do nm=1, nmat - read(*,'(a)') charbuf - charbuf=adjustl(charbuf) - idx=index(charbuf," ") - mtrx(nm)=charbuf(1:idx-1) - rhs(nm)=adjustl(charbuf(idx+1:)) - call psb_bcast(icontxt,mtrx(nm)) - call psb_bcast(icontxt,rhs(nm)) - end do - - else - ! receive parameters - call psb_bcast(icontxt,kmethd) - call psb_bcast(icontxt,eps) - - call psb_bcast(icontxt,afmt) - - call psb_bcast(icontxt,inparms(1:7)) - - ipart = inparms(1) - itmax = inparms(2) - itrace = inparms(3) - istopc = inparms(4) - irst = inparms(5) - irnum = inparms(6) - ntry = inparms(7) - - call psb_bcast(icontxt,nprecs) - allocate(precs(nprecs)) - - do np=1,nprecs - call psb_bcast(icontxt,pdescr) - precs(np)%descr=pdescr - - call psb_bcast(icontxt,lv1) - - call psb_bcast(icontxt,lv2) - - call psb_bcast(icontxt,pparms(1:npparms)) - call psb_bcast(icontxt,omega) - call psb_bcast(icontxt,thr1) - call psb_bcast(icontxt,thr2) - - precs(np)%lv1 = lv1 - precs(np)%lv2 = lv2 - precs(np)%novr = pparms(1) - precs(np)%restr = pparms(2) - precs(np)%prol = pparms(3) - precs(np)%ftype1 = pparms(4) - precs(np)%fill1 = pparms(5) - precs(np)%mltype = pparms(6) - precs(np)%aggr = pparms(7) - precs(np)%smthkind = pparms(8) - precs(np)%cmat = pparms(9) - precs(np)%smthpos = pparms(10) - precs(np)%cslv = pparms(11) - precs(np)%ftype2 = pparms(12) - precs(np)%fill2 = pparms(13) - precs(np)%jswp = pparms(14) - precs(np)%nlev = pparms(15) - precs(np)%omega = omega - precs(np)%thr1 = thr1 - precs(np)%thr2 = thr2 - end do - - - call psb_bcast(icontxt,nmat) - allocate(mtrx(nmat),rhs(nmat)) - - do nm=1,nmat - - call psb_bcast(icontxt,mtrx(nm)) - call psb_bcast(icontxt,rhs(nm)) - - end do - - end if - - end subroutine get_parms - subroutine pr_usage(iout) - integer iout - write(iout, *) ' number of parameters is incorrect!' - write(iout, *) ' use: hb_sample mtrx_file methd prec [ptype & - &itmax istopc itrace]' - write(iout, *) ' where:' - write(iout, *) ' mtrx_file is stored in hb format' - write(iout, *) ' methd may be: cgstab ' - write(iout, *) ' itmax max iterations [500] ' - write(iout, *) ' istopc stopping criterion [1] ' - write(iout, *) ' itrace 0 (no tracing, default) or ' - write(iout, *) ' >= 0 do tracing every itrace' - write(iout, *) ' iterations ' - write(iout, *) ' prec may be: ilu diagsc none' - write(iout, *) ' ptype partition strategy default 0' - write(iout, *) ' 0: block partition ' - end subroutine pr_usage - - -end program zf_bench - - - - - diff --git a/test/pargen/Makefile b/test/pargen/Makefile index bacabfd2..12ab63a7 100644 --- a/test/pargen/Makefile +++ b/test/pargen/Makefile @@ -10,19 +10,21 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBDIR) $(FIFLAG). EXEDIR=./runs -ppde: ppde.o - $(F90LINK) ppde.o -o ppde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) +ppde: ppde.o data_input.o + $(F90LINK) ppde.o data_input.o -o ppde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) /bin/mv ppde $(EXEDIR) -spde: spde.o - $(F90LINK) spde.o -o spde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) +spde: spde.o data_input.o + $(F90LINK) spde.o data_input.o -o spde $(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS) /bin/mv spde $(EXEDIR) +ppde.o spde.o: data_input.o + .f90.o: $(MPF90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $< clean: - /bin/rm -f ppde.o $(EXEDIR)/ppde spde.o $(EXEDIR)/spde + /bin/rm -f data_input.o ppde.o $(EXEDIR)/ppde spde.o $(EXEDIR)/spde verycleanlib: (cd ../..; make veryclean) diff --git a/test/pargen/ppde.f90 b/test/pargen/ppde.f90 index f53d150a..5a74c271 100644 --- a/test/pargen/ppde.f90 +++ b/test/pargen/ppde.f90 @@ -77,51 +77,6 @@ ! ! u(x,y) = rhs(x,y) ! - -module data_input - - interface read_data - module procedure read_char, read_int, read_double - end interface read_data - -contains - - subroutine read_char(val,file) - character(len=*), intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),'(a)') val -!!$ write(0,*) 'read_char got value: "',val,'"' - end subroutine read_char - subroutine read_int(val,file) - integer, intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),*) val -!!$ write(0,*) 'read_int got value: ',val - end subroutine read_int - subroutine read_double(val,file) - use psb_base_mod - real(psb_dpk_), intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),*) val -!!$ write(0,*) 'read_double got value: ',val - end subroutine read_double -end module data_input - program ppde use psb_base_mod use mld_prec_mod @@ -159,7 +114,7 @@ program ppde integer :: novr ! number of overlap layers character(len=16) :: restr ! restriction over application of as character(len=16) :: prol ! prolongation over application of as - character(len=16) :: solve ! Factorization type: ILU, SuperLU, UMFPACK. + character(len=16) :: solve ! Factorization type: ILU, SuperLU, UMFPACK. integer :: fill1 ! Fill-in for factorization 1 real(psb_dpk_) :: thr1 ! Threshold for fact. 1 ILU(T) integer :: nlev ! Number of levels in multilevel prec. @@ -168,11 +123,13 @@ program ppde character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing character(len=16) :: cmat ! coarse mat - character(len=16) :: csolve ! Factorization type: ILU, SuperLU, UMFPACK. + character(len=16) :: csolve ! Coarse solver: bjac, umf, slu, sludist + character(len=16) :: csbsolve ! Coarse subsolver: ILU, ILU(T), SuperLU, UMFPACK. integer :: cfill ! Fill-in for factorization 1 real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_dpk_) :: omega ! smoother omega + real(psb_dpk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prectype ! other variables @@ -234,23 +191,22 @@ program ppde call mld_precset(prec,mld_sub_fillin_,prectype%fill1,info) call mld_precset(prec,mld_fact_thrs_,prectype%thr1,info) if (psb_toupper(prectype%prec) =='ML') then - call mld_precset(prec,mld_aggr_kind_,prectype%aggrkind,info) - call mld_precset(prec,mld_aggr_alg_,prectype%aggr_alg,info) - call mld_precset(prec,mld_ml_type_,prectype%mltype,info) - call mld_precset(prec,mld_ml_type_,prectype%mltype,info) - call mld_precset(prec,mld_smoother_pos_,prectype%smthpos,info) - call mld_precset(prec,mld_coarse_mat_,prectype%cmat,info) - call mld_precset(prec,mld_coarse_solve_,prectype%csolve,info) - call mld_precset(prec,mld_sub_fillin_,prectype%cfill,info,ilev=nlv) - call mld_precset(prec,mld_fact_thrs_,prectype%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smoother_sweeps_,prectype%cjswp,info,ilev=nlv) - call mld_precset(prec,mld_smoother_sweeps_,prectype%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_kind_, prectype%aggrkind,info) + call mld_precset(prec,mld_aggr_alg_, prectype%aggr_alg,info) + call mld_precset(prec,mld_ml_type_, prectype%mltype, info) + call mld_precset(prec,mld_smoother_pos_, prectype%smthpos, info) + call mld_precset(prec,mld_aggr_thresh_, prectype%athres, info) + call mld_precset(prec,mld_coarse_solve_, prectype%csolve, info) + call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) + call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) + call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) + call mld_precset(prec,mld_coarse_fthrs_, prectype%cthres, info) + call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) if (prectype%omega>=0.0) then - call mld_precset(prec,mld_aggr_damp_,prectype%omega,info,ilev=nlv) + call mld_precset(prec,mld_aggr_damp_,prectype%omega,info) end if end if - call psb_barrier(ictxt) t1 = psb_wtime() call mld_precbld(a,desc_a,prec,info) @@ -360,10 +316,12 @@ contains call read_data(prectype%smthpos,5) ! side: pre, post, both smoothing call read_data(prectype%cmat,5) ! coarse mat call read_data(prectype%csolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%csbsolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. call read_data(prectype%cfill,5) ! Fill-in for factorization 1 call read_data(prectype%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prectype%cjswp,5) ! Jacobi sweeps call read_data(prectype%omega,5) ! smoother omega + call read_data(prectype%athres,5) ! smoother aggr thresh end if end if @@ -393,10 +351,12 @@ contains call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat call psb_bcast(ictxt,prectype%csolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%csbsolve) ! Factorization type: ILU, SuperLU, UMFPACK. call psb_bcast(ictxt,prectype%cfill) ! Fill-in for factorization 1 call psb_bcast(ictxt,prectype%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(ictxt,prectype%cjswp) ! Jacobi sweeps call psb_bcast(ictxt,prectype%omega) ! smoother omega + call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh end if if (iam==psb_root_) then diff --git a/test/pargen/spde.f90 b/test/pargen/spde.f90 index 7d26e6c8..110cda67 100644 --- a/test/pargen/spde.f90 +++ b/test/pargen/spde.f90 @@ -77,64 +77,6 @@ ! ! u(x,y) = rhs(x,y) ! - -module data_input - - interface read_data - module procedure read_char, read_int,& - & read_double, read_single - end interface read_data - -contains - - subroutine read_char(val,file) - character(len=*), intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),'(a)') val -!!$ write(0,*) 'read_char got value: "',val,'"' - end subroutine read_char - subroutine read_int(val,file) - integer, intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),*) val -!!$ write(0,*) 'read_int got value: ',val - end subroutine read_int - subroutine read_single(val,file) - use psb_base_mod - real(psb_spk_), intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),*) val -!!$ write(0,*) 'read_double got value: ',val - end subroutine read_single - subroutine read_double(val,file) - use psb_base_mod - real(psb_dpk_), intent(out) :: val - integer, intent(in) :: file - character(len=1024) :: charbuf - integer :: idx - read(file,'(a)')charbuf - charbuf = adjustl(charbuf) - idx=index(charbuf,"!") - read(charbuf(1:idx-1),*) val -!!$ write(0,*) 'read_double got value: ',val - end subroutine read_double -end module data_input - program spde use psb_base_mod use mld_prec_mod @@ -181,11 +123,13 @@ program spde character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing character(len=16) :: cmat ! coarse mat - character(len=16) :: csolve ! Factorization type: ILU, SuperLU, UMFPACK. + character(len=16) :: csolve ! Coarse solver: bjac, umf, slu, sludist + character(len=16) :: csbsolve ! Coarse subsolver: ILU, ILU(T), SuperLU, UMFPACK. integer :: cfill ! Fill-in for factorization 1 real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_spk_) :: omega ! smoother omega + real(psb_spk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prectype ! other variables @@ -247,23 +191,22 @@ program spde call mld_precset(prec,mld_sub_fillin_,prectype%fill1,info) call mld_precset(prec,mld_fact_thrs_,prectype%thr1,info) if (psb_toupper(prectype%prec) =='ML') then - call mld_precset(prec,mld_aggr_kind_,prectype%aggrkind,info) - call mld_precset(prec,mld_aggr_alg_,prectype%aggr_alg,info) - call mld_precset(prec,mld_ml_type_,prectype%mltype,info) - call mld_precset(prec,mld_ml_type_,prectype%mltype,info) - call mld_precset(prec,mld_smoother_pos_,prectype%smthpos,info) - call mld_precset(prec,mld_coarse_mat_,prectype%cmat,info) - call mld_precset(prec,mld_coarse_solve_,prectype%csolve,info) - call mld_precset(prec,mld_sub_fillin_,prectype%cfill,info,ilev=nlv) - call mld_precset(prec,mld_fact_thrs_,prectype%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smoother_sweeps_,prectype%cjswp,info,ilev=nlv) - call mld_precset(prec,mld_smoother_sweeps_,prectype%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_kind_, prectype%aggrkind,info) + call mld_precset(prec,mld_aggr_alg_, prectype%aggr_alg,info) + call mld_precset(prec,mld_ml_type_, prectype%mltype, info) + call mld_precset(prec,mld_smoother_pos_, prectype%smthpos, info) + call mld_precset(prec,mld_aggr_thresh_, prectype%athres, info) + call mld_precset(prec,mld_coarse_solve_, prectype%csolve, info) + call mld_precset(prec,mld_coarse_subsolve_, prectype%csbsolve,info) + call mld_precset(prec,mld_coarse_mat_, prectype%cmat, info) + call mld_precset(prec,mld_coarse_fillin_, prectype%cfill, info) + call mld_precset(prec,mld_coarse_fthrs_, prectype%cthres, info) + call mld_precset(prec,mld_coarse_sweeps_, prectype%cjswp, info) if (prectype%omega>=0.0) then - call mld_precset(prec,mld_aggr_damp_,prectype%omega,info,ilev=nlv) + call mld_precset(prec,mld_aggr_damp_,prectype%omega,info) end if end if - call psb_barrier(ictxt) t1 = psb_wtime() call mld_precbld(a,desc_a,prec,info) @@ -373,10 +316,12 @@ contains call read_data(prectype%smthpos,5) ! side: pre, post, both smoothing call read_data(prectype%cmat,5) ! coarse mat call read_data(prectype%csolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. + call read_data(prectype%csbsolve,5) ! Factorization type: ILU, SuperLU, UMFPACK. call read_data(prectype%cfill,5) ! Fill-in for factorization 1 call read_data(prectype%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prectype%cjswp,5) ! Jacobi sweeps call read_data(prectype%omega,5) ! smoother omega + call read_data(prectype%athres,5) ! smoother aggr thresh end if end if @@ -406,10 +351,12 @@ contains call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat call psb_bcast(ictxt,prectype%csolve) ! Factorization type: ILU, SuperLU, UMFPACK. + call psb_bcast(ictxt,prectype%csbsolve) ! Factorization type: ILU, SuperLU, UMFPACK. call psb_bcast(ictxt,prectype%cfill) ! Fill-in for factorization 1 call psb_bcast(ictxt,prectype%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(ictxt,prectype%cjswp) ! Jacobi sweeps call psb_bcast(ictxt,prectype%omega) ! smoother omega + call psb_bcast(ictxt,prectype%athres) ! smoother aggr thresh end if if (iam==psb_root_) then