fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

mld2p4:

Browse Source 
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.
stopcriterion
Salvatore Filippone 17 years ago
parent 23cdcff940
commit a9722dc049
19 changed files with 2274 additions and 1490 deletions
Show Stats Download Patch File Download Diff File

58
docs/pdf/building.tex
Unescape Escape View File

@ -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}\\

13
docs/pdf/distribution.tex
Unescape Escape View File

@ -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)
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.

4
docs/pdf/license.tex
Unescape Escape View File

@ -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

13
docs/pdf/userinterface.tex
Unescape Escape View File

@ -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|.

40
examples/fileread/Makefile
Unescape Escape View File

@ -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)

91
examples/fileread/data_input.f90
Unescape Escape View File

@ -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

309
examples/fileread/mld_dexample_1lev.f90
Unescape Escape View File

@ -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

348
examples/fileread/mld_dexample_ml.f90
Unescape Escape View File

@ -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

40
examples/pdegen/Makefile
Unescape Escape View File

@ -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)

91
examples/pdegen/data_input.f90
Unescape Escape View File

@ -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

598
examples/pdegen/mld_dexample_1lev.f90
Unescape Escape View File

@ -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<idim) then
val(element)=-b3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x)*idim*idim+(y-1)*idim+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
t3 = psb_wtime()
call psb_spins(element-1,irow,icol,val,a,desc_a,info)
if(info /= 0) exit
tins = tins + (psb_wtime()-t3)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
if(info /= 0) exit
zt(1)=0.d0
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
if(info /= 0) exit
end if
end do
end do
call psb_barrier(ictxt)
t2 = psb_wtime()-t1
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
deallocate(val,irow,icol)
t1 = psb_wtime()
call psb_cdasb(desc_a,info)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_)
call psb_barrier(ictxt)
tasb = psb_wtime()-t1
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
call psb_amx(ictxt,t2)
call psb_amx(ictxt,tins)
call psb_amx(ictxt,tasb)
if(iam == psb_root_) then
write(*,'("The matrix has been generated and assembeld in ",a3," format.")')&
& a%fida(1:3)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2
write(*,'("-assembly time : ",es10.4)')tasb
end if
call psb_geasb(b,desc_a,info)
call psb_geasb(xv,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end subroutine create_matrix
end program mld_dexample_1lev
!
! functions parametrizing the differential equation
!
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a1
real(psb_dpk_) :: x,y,z
!a1=1.d0
a1=0.d0
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a2
real(psb_dpk_) :: x,y,z
!a2=2.d1*y
a2=0.d0
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a3
real(psb_dpk_) :: x,y,z
!a3=1.d0
a3=0.d0
end function a3
function a4(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a4
real(psb_dpk_) :: x,y,z
!a4=1.d0
a4=0.d0
end function a4
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b1
real(psb_dpk_) :: x,y,z
b1=1.d0
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b2
real(psb_dpk_) :: x,y,z
b2=1.d0
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b3
real(psb_dpk_) :: x,y,z
b3=1.d0
end function b3

636
examples/pdegen/mld_dexample_ml.f90
Unescape Escape View File

@ -0,0 +1,636 @@
!!$
!!$
!!$ 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 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_ml
use psb_base_mod
use mld_prec_mod
use psb_krylov_mod
use psb_util_mod
use data_input
implicit none
! input parameters
! sparse matrices
type(psb_dspmat_type) :: A
! sparse matrices descriptor
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 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,amatsize,descsize,precsize
integer :: idim, 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,choice,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(*,'(" ")')
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 solver parameters and 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,choice,idim,itmax,tol)
use psb_base_mod
implicit none
integer :: choice, 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(choice,5)
call read_data(idim,5)
call read_data(itmax,5)
call read_data(tol,5)
end if
call psb_bcast(ictxt,choice)
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<idim) then
val(element)=-b3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x)*idim*idim+(y-1)*idim+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
t3 = psb_wtime()
call psb_spins(element-1,irow,icol,val,a,desc_a,info)
if(info /= 0) exit
tins = tins + (psb_wtime()-t3)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
if(info /= 0) exit
zt(1)=0.d0
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
if(info /= 0) exit
end if
end do
end do
call psb_barrier(ictxt)
t2 = psb_wtime()-t1
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
deallocate(val,irow,icol)
t1 = psb_wtime()
call psb_cdasb(desc_a,info)
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_)
call psb_barrier(ictxt)
tasb = psb_wtime()-t1
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
call psb_amx(ictxt,t2)
call psb_amx(ictxt,tins)
call psb_amx(ictxt,tasb)
if(iam == psb_root_) then
write(*,'("The matrix has been generated and assembeld in ",a3," format.")')&
& a%fida(1:3)
write(*,'("-pspins time : ",es10.4)')tins
write(*,'("-insert time : ",es10.4)')t2
write(*,'("-assembly time : ",es10.4)')tasb
end if
call psb_geasb(b,desc_a,info)
call psb_geasb(xv,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end subroutine create_matrix
end program mld_dexample_ml
!
! functions parametrizing the differential equation
!
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a1
real(psb_dpk_) :: x,y,z
! a1=1.d0
a1=0.d0
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a2
real(psb_dpk_) :: x,y,z
! a2=2.d1*y
a2=0.d0
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a3
real(psb_dpk_) :: x,y,z
! a3=1.d0
a3=0.d0
end function a3
function a4(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a4
real(psb_dpk_) :: x,y,z
! a4=1.d0
a4=0.d0
end function a4
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b1
real(psb_dpk_) :: x,y,z
b1=1.d0
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b2
real(psb_dpk_) :: x,y,z
b2=1.d0
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b3
real(psb_dpk_) :: x,y,z
b3=1.d0
end function b3

14
test/fileread/Makefile
Unescape Escape View File

@ -6,21 +6,14 @@ MLD_LIB=-L$(MLDLIBDIR) -lmld_krylov -lmld_prec
PSBLAS_LIB= -L$(PSBDIR) -lpsb_util -lpsb_base
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBDIR) $(FIFLAG).
DFOBJS=df_bench.o
DFSOBJS=df_sample.o data_input.o
SFSOBJS=sf_sample.o data_input.o
CFSOBJS=cf_sample.o data_input.o
ZFSOBJS=zf_sample.o data_input.o
ZFOBJS=zf_bench.o
EXEDIR=./runs
all: df_bench zf_bench sf_sample df_sample cf_sample zf_sample
df_bench: $(DFOBJS)
$(F90LINK) $(LINKOPT) $(DFOBJS) -o df_bench \
$(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS)
/bin/mv df_bench $(EXEDIR)
all: sf_sample df_sample cf_sample zf_sample
df_sample: $(DFSOBJS)
$(F90LINK) $(LINKOPT) $(DFSOBJS) -o df_sample \
@ -42,11 +35,6 @@ zf_sample: $(ZFSOBJS)
$(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS)
/bin/mv zf_sample $(EXEDIR)
zf_bench: $(ZFOBJS)
$(F90LINK) $(LINKOPT) $(ZFOBJS) -o zf_bench \
$(MLD_LIB) $(PSBLAS_LIB) $(LDLIBS)
/bin/mv zf_bench $(EXEDIR)
sf_sample.o: data_input.o
df_sample.o: data_input.o
cf_sample.o: data_input.o

648
test/fileread/df_bench.f90
Unescape Escape View File

@ -1,648 +0,0 @@
program df_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_dspmat_type) :: a, aux_a
! preconditioner data
type(mld_dprec_type) :: pre
integer :: igsmth, matop, novr
! dense matrices
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
real(psb_dpk_), allocatable , save :: b_col(:), x_col(:), r_col(:), &
& x_col_glob(:), r_col_glob(:)
real(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='df_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,psb_root_)
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),psb_root_)
else
call psb_bcast(ictxt,m_problem,psb_root_)
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),psb_root_)
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))
if (.true.) then
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
call psb_matdist(aux_a, a, part_block, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt)
end if
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(0,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
!!$ open(20+iam)
!!$ call psb_cdprt(20+iam,desc_a,short=.false.)
!!$ close(20+iam)
!!$ write(0,*) iam,'After CDPRT '
!!$ call flush(0)
!!$ call flush(6)
!!$ call psb_barrier(ictxt)
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(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

84
test/fileread/runs/drt.sh
Unescape Escape View File

@ -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 <<EOF
mpirun -np $np -machinefile locm df_bench >>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <<EOF
out.${np}p.$date Out file 1: summary
stat.${np}p.$date Out file 2: detailed for statistics
BICGSTAB iterative method to use
1.D-6 EPS
CSR Matrix format
0 IPART: Partition method: 0: BLOCK 1: BLK2 2:GRAPH
00250 ITMAX
-1 ITRACE
2 ISTOPC 1: NBE Infty 2: |r|2/|b|2
30 IRST Restart parameter for GMRES and BiCGSTAB(L)
0 RENUM: 0: none 1: global indices (2: GPS band reduction)
$ntry NTRY for each comb. print out best timings
30 NPRCS nov rst prl fc1 fl1 mlt agg smt cm smp csv ft2 fl2 jsw nl omg th1 th2 name
none none 0 0 0 0 0 0 0 0 0 0 2 0 0 0 1 -1.0 1e-4 1e-4 NOPREC
diag none 0 0 0 1 0 2 0 1 0 2 2 1 0 4 1 -1.0 1e-4 1e-4 DIAG
bjac none 0 0 0 1 0 2 0 1 0 2 2 1 0 4 1 -1.0 1e-4 1e-4 BJAC
as none 0 1 0 1 0 2 0 1 0 2 2 1 0 4 1 -1.0 1e-4 1e-4 RAS
as none 1 1 0 1 0 2 0 1 0 2 2 1 0 4 1 -1.0 1e-4 1e-4 RAS
as none 2 1 0 1 0 2 0 1 0 2 2 1 0 4 1 -1.0 1e-4 1e-4 RAS
as ml 0 1 0 1 0 2 0 1 0 2 2 1 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-I-D4
as ml 1 1 0 1 0 2 0 1 0 2 2 1 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-I-D4
as ml 2 1 0 1 0 2 0 1 0 2 2 1 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-I-D4
as ml 0 1 0 1 0 2 0 1 0 2 2 5 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-U-D4
as ml 1 1 0 1 0 2 0 1 0 2 2 5 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-U-D4
as ml 2 1 0 1 0 2 0 1 0 2 2 5 0 4 2 -1.0 1e-4 1e-4 2L-M-RAS-U-D4
as ml 0 1 0 1 0 2 0 1 0 2 2 1 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-I-D4
as ml 1 1 0 1 0 2 0 1 0 2 2 1 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-I-D4
as ml 2 1 0 1 0 2 0 1 0 2 2 1 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-I-D4
as ml 0 1 0 1 0 2 0 1 0 2 2 5 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-U-D4
as ml 1 1 0 1 0 2 0 1 0 2 2 5 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-U-D4
as ml 2 1 0 1 0 2 0 1 0 2 2 5 0 4 3 -1.0 1e-4 1e-4 3L-M-RAS-U-D4
as ml 0 1 0 1 0 2 0 1 1 2 2 1 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-I-R
as ml 1 1 0 1 0 2 0 1 1 2 2 1 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-I-R
as ml 2 1 0 1 0 2 0 1 1 2 2 1 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-I-R
as ml 0 1 0 1 0 2 0 1 1 2 2 5 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-U-R
as ml 1 1 0 1 0 2 0 1 1 2 2 5 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-U-R
as ml 2 1 0 1 0 2 0 1 1 2 2 5 0 1 2 -1.0 1e-4 1e-4 2L-M-RAS-U-R
as ml 0 1 0 1 0 2 0 1 1 2 2 1 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-I-R
as ml 1 1 0 1 0 2 0 1 1 2 2 1 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-I-R
as ml 2 1 0 1 0 2 0 1 1 2 2 1 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-I-R
as ml 0 1 0 1 0 2 0 1 1 2 2 5 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-U-R
as ml 1 1 0 1 0 2 0 1 1 2 2 5 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-U-R
as ml 2 1 0 1 0 2 0 1 1 2 2 5 0 1 3 -1.0 1e-4 1e-4 3L-M-RAS-U-R
2 Number of matrices
kivap004.mtx none
kivap001.mtx none
thm50x30.mtx none
thm200x120.mtx none
thm1000x600.mtx none
a400x400.mtx b400.mtx
kivap007.mtx none
!!! preconditioner templates
bja none 0 0 0 0 0 0 0 0 0 0 0 0.0 Block Jacobi
none none 0 0 0 0 0 0 0 0 0 0 0 0.1 No preconditioner
diagsc none 0 0 0 0 0 0 0 0 0 0 0 0.0 Diagonal scaling
as none 1 4 1 0 0 0 0 0 0 0 0 0.0 Additive Schwarz 1 overlap
as none 1 4 0 0 0 0 0 0 0 0 0 0.0 Restricted Additive Schwarz 1 overlap
EOF
cat out.${np}p.$date >>dat.out.kiva.$date
cat stat.${np}p.$date >>dat.stat.kiva.$date
done

594
test/fileread/zf_bench.f90
Unescape Escape View File

@ -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

12
test/pargen/Makefile
Unescape Escape View File

@ -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)

80
test/pargen/ppde.f90
Unescape Escape View File

@ -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

91
test/pargen/spde.f90
Unescape Escape View File

@ -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

Write Preview
Loading…
Cancel
Save