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

mld2p4-2:

Browse Source 
Makefile
 config/pac.m4
 configure.ac
 configure
 mlprec/Makefile
 mlprec/mld_daggrmap_bld.f90
 mlprec/mld_daggrmat_asb.f90
 mlprec/mld_daggrmat_minnrg_asb.F90
 mlprec/mld_daggrmat_nosmth_asb.F90
 mlprec/mld_daggrmat_smth_asb.F90
 mlprec/mld_das_aply.f90
 mlprec/mld_das_bld.f90
 mlprec/mld_dbaseprec_bld.f90
 mlprec/mld_dcoarse_bld.f90
 mlprec/mld_ddiag_bld.f90
 mlprec/mld_dfact_bld.f90
 mlprec/mld_dilu0_fact.f90
 mlprec/mld_dilu_bld.f90
 mlprec/mld_diluk_fact.f90
 mlprec/mld_dilut_fact.f90
 mlprec/mld_dmlprec_bld.f90
 mlprec/mld_dprecaply.f90
 mlprec/mld_dprecbld.f90
 mlprec/mld_dslu_bld.f90
 mlprec/mld_dslud_bld.f90
 mlprec/mld_dsp_renum.f90
 mlprec/mld_dsub_solve.f90
 mlprec/mld_dumf_bld.f90
 mlprec/mld_inner_mod.f90
 mlprec/mld_prec_mod.f90
 mlprec/mld_prec_type.f90
 tests/pdegen/Makefile
 tests/pdegen/ppde.f90

MLD2P4: first compilable version for D, still a lot to do to make it
RUN properly.
stopcriterion
Salvatore Filippone 15 years ago
parent 5dfc73dcad
commit a64fad80c6
33 changed files with 4882 additions and 2976 deletions
Show Stats Download Patch File Download Diff File

2
Makefile
Unescape Escape View File

@ -2,7 +2,7 @@ include Make.inc
all: library
library: libdir mlp kryl
library: libdir mlp
libdir:
(if test ! -d lib ; then mkdir lib; fi)

49
config/pac.m4
Unescape Escape View File

@ -694,20 +694,65 @@ if test "x$mld2p4_cv_umfpackdir" != "x"; then
LIBS="-L$mld2p4_cv_umfpackdir $LIBS"
UMF_INCLUDES="-I$mld2p4_cv_umfpackdir"
CPPFLAGS="$UMF_INCLUDES $CPPFLAGS"
UMF_LIBS="-L$mld2p4_cv_umfpackdir"
UMF_LIBDIR="-L$mld2p4_cv_umfpackdir"
fi
AC_MSG_NOTICE([umfp dir $mld2p4_cv_umfpackdir])
AC_CHECK_HEADER([umfpack.h],
[pac_umf_header_ok=yes],
[pac_umf_header_ok=no; UMF_INCLUDES=""])
if test "x$pac_umf_header_ok" == "xno" ; then
dnl Maybe Include or include subdirs?
unset ac_cv_header_umfpack_h
UMF_INCLUDES="-I$mld2p4_cv_umfpackdir/include -I$mld2p4_cv_umfpackdir/Include "
CPPFLAGS="$UMF_INCLUDES $SAVE_CPPFLAGS"
AC_MSG_CHECKING([for umfpack_di_symbolic in $UMF_INCLUDES])
AC_CHECK_HEADER([umfpack.h],
[pac_umf_header_ok=yes],
[pac_umf_header_ok=no; UMF_INCLUDES=""])
fi
if test "x$pac_umf_header_ok" == "xno" ; then
dnl Maybe new structure with UMFPACK UFconfig AMD?
unset ac_cv_header_umfpack_h
UMF_INCLUDES="-I$mld2p4_cv_umfpackdir/UFconfig -I$mld2p4_cv_umfpackdir/UMFPACK/Include -I$mld2p4_cv_umfpackdir/AMD/Include"
CPPFLAGS="$UMF_INCLUDES $SAVE_CPPFLAGS"
AC_CHECK_HEADER([umfpack.h],
[pac_umf_header_ok=yes],
[pac_umf_header_ok=no; UMF_INCLUDES=""])
fi
if test "x$pac_umf_header_ok" == "xyes" ; then
UMF_LIBS="$mld2p4_cv_umfpack $UMF_LIBS"
UMF_LIBS="$mld2p4_cv_umfpack $UMF_LIBDIR"
LIBS="$UMF_LIBS -lm $LIBS";
AC_MSG_CHECKING([for umfpack_di_symbolic in $UMF_LIBS])
AC_TRY_LINK_FUNC(umfpack_di_symbolic,
[mld2p4_cv_have_umfpack=yes;pac_umf_lib_ok=yes; ],
[mld2p4_cv_have_umfpack=no;pac_umf_lib_ok=no; UMF_LIBS=""])
AC_MSG_RESULT($pac_umf_lib_ok)
if test "x$pac_umf_lib_ok" == "xno" ; then
dnl Maybe Lib or lib?
UMF_LIBDIR="-L$mld2p4_cv_umfpackdir/Lib -L$mld2p4_cv_umfpackdir/lib"
UMF_LIBS="$mld2p4_cv_umfpack $UMF_LIBDIR -lm $SAVE_LIBS"
LIBS="$UMF_LIBS"
AC_MSG_CHECKING([for umfpack_di_symbolic in $UMF_LIBS])
AC_TRY_LINK_FUNC(umfpack_di_symbolic,
[mld2p4_cv_have_umfpack=yes;pac_umf_lib_ok=yes; ],
[mld2p4_cv_have_umfpack=no;pac_umf_lib_ok=no; UMF_LIBS=""])
AC_MSG_RESULT($pac_umf_lib_ok)
fi
if test "x$pac_umf_lib_ok" == "xno" ; then
dnl Maybe UMFPACK/Lib?
UMF_LIBDIR="-L$mld2p4_cv_umfpackdir/AMD/Lib -L$mld2p4_cv_umfpackdir/UMFPACK/Lib"
UMF_LIBS="$mld2p4_cv_umfpack $UMF_LIBDIR -lm $SAVE_LIBS"
LIBS="$UMF_LIBS"
AC_MSG_CHECKING([for umfpack_di_symbolic in $UMF_LIBS])
AC_TRY_LINK_FUNC(umfpack_di_symbolic,
[mld2p4_cv_have_umfpack=yes;pac_umf_lib_ok=yes; ],
[mld2p4_cv_have_umfpack=no;pac_umf_lib_ok=no; UMF_LIBS=""])
AC_MSG_RESULT($pac_umf_lib_ok)
fi
fi
LIBS="$SAVE_LIBS";
CPPFLAGS="$SAVE_CPPFLAGS";

2143
configure vendored
View File

File diff suppressed because it is too large Load Diff

46
configure.ac
Unescape Escape View File

@ -38,11 +38,11 @@ dnl NOTE : odd configurations like ifc + gcc still await in the mist of the unkn
###############################################################################
# NOTE: the literal for version (the second argument to AC_INIT should be a literal!)
AC_INIT([MLD2P4],1.1, bugreport@mld2p4.it)
AC_INIT([MLD2P4],2.0, bugreport@mld2p4.it)
# VERSION is the file containing the PSBLAS version code
# FIXME
mld2p4_cv_version="1.1"
mld2p4_cv_version="2.0"
# A sample source file
AC_CONFIG_SRCDIR([mlprec/mld_prec_type.f90])
@ -116,7 +116,7 @@ AC_MSG_RESULT([$INSTALL_DIR $INSTALL_INCLUDEDIR $INSTALL_LIBDIR $INSTALL_DOCSDIR
# Compilers detection: FC,F77,CC should be set, if found.
###############################################################################
AC_PROG_FC([xlf95 xlf90 xlf pgf95 pgf90 ifort ifc gfortran])
AC_PROG_FC([xlf95 xlf90 xlf pgf95 pgf90 ifort ifc nagfor gfortran])
AC_PROG_CC([xlc pgcc icc gcc ])
dnl AC_PROG_CXX
@ -136,12 +136,13 @@ if eval "$FC -qversion 2>&1 | grep XL 2>/dev/null" ; then
# More problems could be undocumented yet.
fi
PAC_ARG_WITH_LIBS
PAC_ARG_WITH_FLAGS(clibs,CLIBS)
PAC_ARG_WITH_FLAGS(flibs,FLIBS)
PAC_ARG_WITH_FLAGS(library-path,LIBRARYPATH)
PAC_ARG_WITH_FLAGS(include-path,INCLUDEPATH)
PAC_ARG_WITH_FLAGS(module-path,MODULE_PATH)
dnl PAC_ARG_WITH_FLAGS(clibs,CLIBS)
dnl PAC_ARG_WITH_FLAGS(flibs,FLIBS)
dnl PAC_ARG_WITH_FLAGS(library-path,LIBRARYPATH)
dnl PAC_ARG_WITH_FLAGS(include-path,INCLUDEPATH)
dnl PAC_ARG_WITH_FLAGS(module-path,MODULE_PATH)
dnl CPPFLAGS="$INCLUDEPATH $CPPFLAGS"
###############################################################################
# BLAS library presence checks
###############################################################################
@ -176,12 +177,41 @@ else
UMF_FLAGS=""
fi
PAC_ARG_SERIAL_MPI
#Note : we miss the name of the Intel C compiler
if test x"$pac_cv_serial_mpi" == x"yes" ; then
FAKEMPI="fakempi.o";
MPIFC="$FC";
MPIF77="$F77";
MPICC="$CC";
else
AC_LANG([C])
if test "X$MPICC" = "X" ; then
# This is our MPICC compiler preference: it will override ACX_MPI's first try.
AC_CHECK_PROGS([MPICC],[mpxlc mpcc mpicc])
fi
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for C]])])
AC_LANG(Fortran 77)
if test "X$MPIF77" = "X" ; then
# This is our MPIFC compiler preference: it will override ACX_MPI's first try.
AC_CHECK_PROGS([MPIF77],[mpxlf mpf77 mpif77])
fi
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for Fortran 77]])])
AC_LANG([Fortran])
if test "X$MPIFC" = "X" ; then
# This is our MPIFC compiler preference: it will override ACX_MPI's first try.
AC_CHECK_PROGS([MPIFC],[mpxlf95 mpxlf90 mpf95 mpf90 mpif95 mpif90 ])
fi
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for Fortran]])])
fi
# We leave a default language for the next checks.
dnl AC_LANG([Fortran 77])
AC_LANG([C])
PAC_CHECK_SUPERLU
if test "x$mld2p4_cv_have_superlu" == "xyes" ; then

92
mlprec/Makefile
Unescape Escape View File

@ -7,49 +7,65 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBLIBDIR)
MODOBJS=mld_prec_type.o mld_prec_mod.o mld_inner_mod.o mld_move_alloc_mod.o
MPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o \
mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o\
mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o \
mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o
MPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o
MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o
INNEROBJS=mld_scoarse_bld.o mld_dcoarse_bld.o \
mld_ccoarse_bld.o mld_zcoarse_bld.o \
mld_smlprec_bld.o mld_dmlprec_bld.o mld_cmlprec_bld.o mld_zmlprec_bld.o\
mld_sas_bld.o mld_sslu_bld.o mld_sumf_bld.o mld_silu0_fact.o\
mld_ssp_renum.o mld_sfact_bld.o mld_silu_bld.o \
mld_sbaseprec_bld.o mld_sdiag_bld.o mld_saggrmap_bld.o \
mld_smlprec_aply.o mld_sslud_bld.o\
mld_sbaseprec_aply.o mld_ssub_aply.o mld_ssub_solve.o \
mld_sas_aply.o mld_saggrmat_asb.o \
mld_das_bld.o mld_dslu_bld.o mld_dumf_bld.o mld_dilu0_fact.o\
mld_dsp_renum.o mld_dfact_bld.o mld_dilu_bld.o \
INNEROBJS= mld_dcoarse_bld.o \
mld_dmlprec_bld.o\
mld_das_bld.o mld_dslu_bld.o mld_dumf_bld.o \
mld_dilu_bld.o mld_dilu0_fact.o mld_dfact_bld.o \
mld_diluk_fact.o mld_dilut_fact.o \
mld_dbaseprec_bld.o mld_ddiag_bld.o mld_daggrmap_bld.o \
mld_dmlprec_aply.o mld_dslud_bld.o\
mld_dbaseprec_aply.o mld_dsub_aply.o mld_dsub_solve.o \
mld_das_aply.o mld_daggrmat_asb.o \
mld_cas_bld.o mld_cslu_bld.o mld_cumf_bld.o mld_cilu0_fact.o\
mld_csp_renum.o mld_cfact_bld.o mld_cilu_bld.o \
mld_cbaseprec_bld.o mld_cdiag_bld.o mld_caggrmap_bld.o \
mld_cmlprec_aply.o mld_cslud_bld.o\
mld_cbaseprec_aply.o mld_csub_aply.o mld_csub_solve.o \
mld_cas_aply.o mld_caggrmat_asb.o\
mld_zas_bld.o mld_zslu_bld.o mld_zumf_bld.o mld_zilu0_fact.o\
mld_zsp_renum.o mld_zfact_bld.o mld_zilu_bld.o \
mld_zbaseprec_bld.o mld_zdiag_bld.o mld_zaggrmap_bld.o \
mld_zmlprec_aply.o mld_zslud_bld.o\
mld_zbaseprec_aply.o mld_zsub_aply.o mld_zsub_solve.o \
mld_zas_aply.o mld_zaggrmat_asb.o\
mld_siluk_fact.o mld_ciluk_fact.o mld_silut_fact.o mld_cilut_fact.o \
mld_diluk_fact.o mld_ziluk_fact.o mld_dilut_fact.o mld_zilut_fact.o \
mld_das_aply.o mld_daggrmat_asb.o mld_dsp_renum.o \
$(MPFOBJS)
OUTEROBJS=mld_sprecbld.o mld_sprecset.o mld_sprecinit.o\
mld_sprecaply.o \
mld_dprecbld.o mld_dprecset.o mld_dprecinit.o\
mld_dprecaply.o \
mld_cprecbld.o mld_cprecset.o mld_cprecinit.o\
mld_cprecaply.o \
mld_zprecbld.o mld_zprecset.o mld_zprecinit.o \
mld_zprecaply.o
#
# MPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o \
# mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o\
# mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o \
# mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o
# MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zslud_interface.o
# INNEROBJS=mld_scoarse_bld.o mld_dcoarse_bld.o \
# mld_ccoarse_bld.o mld_zcoarse_bld.o \
# mld_smlprec_bld.o mld_dmlprec_bld.o mld_cmlprec_bld.o mld_zmlprec_bld.o\
# mld_sas_bld.o mld_sslu_bld.o mld_sumf_bld.o mld_silu0_fact.o\
# mld_ssp_renum.o mld_sfact_bld.o mld_silu_bld.o \
# mld_sbaseprec_bld.o mld_sdiag_bld.o mld_saggrmap_bld.o \
# mld_smlprec_aply.o mld_sslud_bld.o\
# mld_sbaseprec_aply.o mld_ssub_aply.o mld_ssub_solve.o \
# mld_sas_aply.o mld_saggrmat_asb.o \
# mld_das_bld.o mld_dslu_bld.o mld_dumf_bld.o mld_dilu0_fact.o\
# mld_dsp_renum.o mld_dfact_bld.o mld_dilu_bld.o \
# mld_dbaseprec_bld.o mld_ddiag_bld.o mld_daggrmap_bld.o \
# mld_dmlprec_aply.o mld_dslud_bld.o\
# mld_dbaseprec_aply.o mld_dsub_aply.o mld_dsub_solve.o \
# mld_das_aply.o mld_daggrmat_asb.o \
# mld_cas_bld.o mld_cslu_bld.o mld_cumf_bld.o mld_cilu0_fact.o\
# mld_csp_renum.o mld_cfact_bld.o mld_cilu_bld.o \
# mld_cbaseprec_bld.o mld_cdiag_bld.o mld_caggrmap_bld.o \
# mld_cmlprec_aply.o mld_cslud_bld.o\
# mld_cbaseprec_aply.o mld_csub_aply.o mld_csub_solve.o \
# mld_cas_aply.o mld_caggrmat_asb.o\
# mld_zas_bld.o mld_zslu_bld.o mld_zumf_bld.o mld_zilu0_fact.o\
# mld_zsp_renum.o mld_zfact_bld.o mld_zilu_bld.o \
# mld_zbaseprec_bld.o mld_zdiag_bld.o mld_zaggrmap_bld.o \
# mld_zmlprec_aply.o mld_zslud_bld.o\
# mld_zbaseprec_aply.o mld_zsub_aply.o mld_zsub_solve.o \
# mld_zas_aply.o mld_zaggrmat_asb.o\
# mld_siluk_fact.o mld_ciluk_fact.o mld_silut_fact.o mld_cilut_fact.o \
# mld_diluk_fact.o mld_ziluk_fact.o mld_dilut_fact.o mld_zilut_fact.o \
# $(MPFOBJS)
OUTEROBJS=mld_dprecbld.o mld_dprecset.o mld_dprecinit.o\
mld_dprecaply.o
#OUTEROBJS=mld_sprecbld.o mld_sprecset.o mld_sprecinit.o\
# mld_sprecaply.o \
# mld_dprecbld.o mld_dprecset.o mld_dprecinit.o\
# mld_dprecaply.o \
# mld_cprecbld.o mld_cprecset.o mld_cprecinit.o\
# mld_cprecaply.o \
# mld_zprecbld.o mld_zprecset.o mld_zprecinit.o \
# mld_zprecaply.o
F90OBJS=$(OUTEROBJS) $(INNEROBJS)
COBJS= mld_sslu_interface.o mld_sumf_interface.o \

62
mlprec/mld_daggrmap_bld.f90
Unescape Escape View File

@ -87,17 +87,17 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
implicit none
! Arguments
integer, intent(in) :: aggr_type
real(psb_dpk_), intent(in) :: theta
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
integer, intent(in) :: aggr_type
real(psb_dpk_), intent(in) :: theta
type(psb_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:)
integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m
type(psb_dspmat_type) :: atmp, atrans
type(psb_d_sparse_mat) :: atmp, atrans
logical :: recovery
integer :: debug_level, debug_unit
integer :: ictxt,np,me,err_act
@ -121,19 +121,20 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info)
call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info)
case (mld_sym_dec_aggr_)
nr = psb_sp_get_nrows(a)
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
nr = a%get_nrows()
call a%csclip(atmp,info,imax=nr,jmax=nr,&
& rscale=.false.,cscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_transp(atmp,atrans,fmt='COO')
call atmp%set_nrows(nr)
call atmp%set_ncols(nr)
if (info == 0) call atrans%transp(atmp)
if (info == 0) call atrans%cscnv(info,type='COO')
if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
atmp%m=nr
atmp%k=nr
if (info == 0) call psb_sp_free(atrans,info)
if (info == 0) call psb_spcnv(atmp,info,afmt='csr')
call atmp%set_nrows(nr)
call atmp%set_ncols(nr)
if (info == 0) call atrans%free()
if (info == 0) call atmp%cscnv(info,type='CSR')
if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info)
if (info == 0) call psb_sp_free(atmp,info)
if (info == 0) call atmp%free()
case default
@ -169,11 +170,11 @@ contains
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
type(psb_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
real(psb_dpk_), intent(in) :: theta
integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:)
integer, intent(out) :: info
! Local variables
integer, allocatable :: ils(:), neigh(:), irow(:), icol(:)
@ -186,7 +187,7 @@ contains
integer :: nrow, ncol, n_ne
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
if (psb_get_errstatus() /= 0) return
info=0
name = 'mld_dec_map_bld'
call psb_erractionsave(err_act)
@ -198,7 +199,7 @@ contains
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
nr = a%m
nr = a%get_nrows()
allocate(ilaggr(nr),neigh(nr),stat=info)
if(info /= 0) then
info=4025
@ -214,7 +215,7 @@ contains
& a_err='real(psb_dpk_)')
goto 9999
end if
call psb_sp_getdiag(a,diag,info)
call a%get_diag(diag,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getdiag')
@ -247,10 +248,10 @@ contains
naggr = naggr + 1
ilaggr(i) = naggr
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
call a%csget(i,i,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getrow')
call psb_errpush(info,name,a_err='csget')
goto 9999
end if
@ -268,7 +269,7 @@ contains
!
! 2. Untouched neighbours of these nodes are marked <0.
!
call psb_neigh(a,i,neigh,n_ne,info,lev=2)
call a%get_neigh(i,neigh,n_ne,info,lev=2)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_neigh')
@ -288,8 +289,7 @@ contains
enddo
if (debug_level >= psb_debug_outer_) then
write(debug_unit,*) me,' ',trim(name),&
& ' Check 1:',count(ilaggr == -(nr+1)),&
& (a%ia1(i),i=a%ia2(1),a%ia2(2)-1)
& ' Check 1:',count(ilaggr == -(nr+1))
end if
!
@ -336,7 +336,7 @@ contains
isz = nr+1
ia = -1
cpling = dzero
call psb_sp_getrow(i,a,nz,irow,icol,val,info)
call a%csget(i,i,nz,irow,icol,val,info)
if (info/=0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getrow')

10
mlprec/mld_daggrmat_asb.f90
Unescape Escape View File

@ -106,11 +106,11 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_donelev_type), intent(inout), target :: p
integer, intent(out) :: info
type(psb_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_donelev_type), intent(inout), target :: p
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me, err_act, icomm

1516
mlprec/mld_daggrmat_minnrg_asb.F90
View File

File diff suppressed because it is too large Load Diff

131
mlprec/mld_daggrmat_nosmth_asb.F90
Unescape Escape View File

@ -61,7 +61,7 @@
!
!
! Arguments:
! a - type(psb_dspmat_type), input.
! a - type(psb_d_sparse_mat), input.
! The sparse matrix structure containing the local part of
! the fine-level matrix.
! desc_a - type(psb_desc_type), input.
@ -93,19 +93,20 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
include 'mpif.h'
#endif
! Arguments
type(psb_dspmat_type), intent(in) :: a
! Arguments
type(psb_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_donelev_type), intent(inout), target :: p
integer, intent(out) :: info
! Local variables
! Local variables
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type) :: b
type(psb_d_sparse_mat) :: b
integer, allocatable :: nzbr(:), idisp(:)
type(psb_dspmat_type) :: am1,am2
type(psb_d_sparse_mat) :: am1,am2
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzt, naggrm1, i
@ -114,7 +115,6 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
info=0
call psb_erractionsave(err_act)
call psb_nullify_sp(b)
ictxt = psb_cd_get_context(desc_a)
icomm = psb_cd_get_mpic(desc_a)
@ -123,9 +123,6 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
@ -152,47 +149,33 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
call psb_sp_all(ncol,ntaggr,am1,ncol,info)
call acoo1%allocate(ncol,ntaggr,ncol)
else
call psb_sp_all(ncol,naggr,am1,ncol,info)
end if
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
call acoo1%allocate(ncol,naggr,ncol)
end if
do i=1,nrow
am1%aspk(i) = done
am1%ia1(i) = i
am1%ia2(i) = ilaggr(i)
acoo1%val(i) = done
acoo1%ia(i) = i
acoo1%ja(i) = ilaggr(i)
end do
am1%infoa(psb_nnz_) = nrow
call psb_spcnv(am1,info,afmt='csr',dupl=psb_dupl_add_)
call psb_transp(am1,am2)
call acoo1%set_nzeros(nrow)
call acoo2%transp(acoo1)
call a%csclip(bcoo,info,jmax=nrow)
call psb_sp_clip(a,b,info,jmax=nrow)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclip')
goto 9999
end if
! Out from sp_clip is always in COO, but just in case..
if (psb_tolower(b%fida) /= 'coo') then
call psb_errpush(4010,name,a_err='spclip NOT COO')
goto 9999
end if
nzt = psb_sp_get_nnzeros(b)
nzt = bcoo%get_nzeros()
do i=1, nzt
b%ia1(i) = ilaggr(b%ia1(i))
b%ia2(i) = ilaggr(b%ia2(i))
bcoo%ia(i) = ilaggr(bcoo%ia(i))
bcoo%ja(i) = ilaggr(bcoo%ja(i))
enddo
b%m = naggr
b%k = naggr
! This is to minimize data exchange
call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
call bcoo%set_nrows(naggr)
call bcoo%set_ncols(naggr)
call bcoo%set_dupl(psb_dupl_add_)
call bcoo%fix(info)
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
@ -206,81 +189,73 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nzbr(me+1) = nzt
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,p%ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_all')
goto 9999
end if
call ac_coo%allocate(ntaggr,ntaggr,nzac)
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,p%ac%aspk,nzbr,idisp,&
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(b%ia1,ndx,mpi_integer,p%ac%ia1,nzbr,idisp,&
call mpi_allgatherv(bcoo%ia,ndx,mpi_integer,ac_coo%ia,nzbr,idisp,&
& mpi_integer,icomm,info)
call mpi_allgatherv(b%ia2,ndx,mpi_integer,p%ac%ia2,nzbr,idisp,&
call mpi_allgatherv(bcoo%ja,ndx,mpi_integer,ac_coo%ja,nzbr,idisp,&
& mpi_integer,icomm,info)
if(info /= 0) then
info=-1
call psb_errpush(info,name)
goto 9999
end if
p%ac%m = ntaggr
p%ac%k = ntaggr
p%ac%infoa(psb_nnz_) = nzac
p%ac%fida='COO'
p%ac%descra='GUN'
call psb_spcnv(p%ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
goto 9999
end if
call ac_coo%set_nzeros(nzac)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call p%ac%mv_from(ac_coo)
else if (p%iprcparm(mld_coarse_mat_) == mld_distr_mat_) then
call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(p%desc_ac,info)
if (info == 0) call psb_sp_clone(b,p%ac,info)
call p%ac%mv_from(bcoo)
if(info /= 0) then
call psb_errpush(4001,name,a_err='Build ac, desc_ac')
goto 9999
end if
call psb_sp_free(b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
else
info = 4001
call psb_errpush(4001,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
else
info = 4001
call psb_errpush(4001,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
call bcoo%free()
deallocate(nzbr,idisp)
call psb_spcnv(p%ac,info,afmt='csr',dupl=psb_dupl_add_)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
call psb_errpush(4010,name,a_err='cscnv')
goto 9999
end if
call am1%mv_from(acoo1)
call am1%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == 0) call am2%mv_from(acoo2)
if (info == 0) call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
!
! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
if (info == 0) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr)
if (info == 0) call psb_sp_free(am1,info)
if (info == 0) call psb_sp_free(am2,info)
if (info == 0) call am1%free()
if (info == 0) call am2%free()
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_Free')
call psb_errpush(4010,name,a_err='linmap build')
goto 9999
end if

432
mlprec/mld_daggrmat_smth_asb.F90
Unescape Escape View File

@ -111,21 +111,22 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
#endif
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_d_sparse_mat), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_donelev_type), intent(inout), target :: p
integer, intent(out) :: info
! Local variables
type(psb_dspmat_type) :: b
type(psb_d_sparse_mat) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, af
type(psb_dspmat_type) :: am3,am4
type(psb_d_sparse_mat) :: am1,am2, am3, am4
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, acoof, acoo3,acoo4, bcoo, cootmp
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsrf, acsr3,acsr4, bcsr
real(psb_dpk_), allocatable :: adiag(:)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
@ -145,14 +146,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_info(ictxt, me, np)
call psb_nullify_sp(b)
call psb_nullify_sp(am3)
call psb_nullify_sp(am4)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
call psb_nullify_sp(AF)
nglob = psb_cd_get_global_rows(desc_a)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
@ -201,7 +194,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
! Get the diagonal D
call psb_sp_getdiag(a,adiag,info)
call a%get_diag(adiag,info)
if (info == 0) &
& call psb_halo(adiag,desc_a,info)
@ -211,44 +204,28 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
! 1. Allocate Ptilde in sparse matrix form
am4%fida='COO'
am4%m=ncol
if (ml_global_nmb) then
am4%k=ntaggr
call psb_sp_all(ncol,ntaggr,am4,ncol,info)
else
am4%k=naggr
call psb_sp_all(ncol,naggr,am4,ncol,info)
endif
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
end if
if (ml_global_nmb) then
call acoo4%allocate(ncol,ntaggr,ncol)
do i=1,ncol
am4%aspk(i) = done
am4%ia1(i) = i
am4%ia2(i) = ilaggr(i)
acoo4%val(i) = done
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
end do
am4%infoa(psb_nnz_) = ncol
else
call acoo4%set_nzeros(ncol)
else
call acoo4%allocate(ncol,naggr,ncol)
do i=1,nrow
am4%aspk(i) = done
am4%ia1(i) = i
am4%ia2(i) = ilaggr(i)
acoo4%val(i) = done
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
end do
am4%infoa(psb_nnz_) = nrow
call acoo4%set_nzeros(nrow)
endif
call acoo4%set_dupl(psb_dupl_add_)
call acsr4%mv_from_coo(acoo4,info)
if (info==0) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
call psb_spcnv(am4,info,afmt='csr',dupl=psb_dupl_add_)
if (info==0) call psb_spcnv(a,am3,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies done.'
@ -257,39 +234,38 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Build the filtered matrix Af from A
!
call psb_spcnv(a,af,info,afmt='csr',dupl=psb_dupl_add_)
if (info==0) call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = dzero
jd = -1
do j=af%ia2(i),af%ia2(i+1)-1
if (af%ia1(j) == i) jd = j
if (abs(af%aspk(j)) < theta*sqrt(abs(adiag(i)*adiag(af%ia1(j))))) then
tmp=tmp+af%aspk(j)
af%aspk(j)=dzero
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) jd = j
if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
tmp=tmp+acsrf%val(j)
acsrf%val(j)=dzero
endif
enddo
if (jd == -1) then
write(0,*) 'Wrong input: we need the diagonal!!!!', i
else
af%aspk(jd)=af%aspk(jd)-tmp
acsrf%val(jd)=acsrf%val(jd)-tmp
end if
enddo
! Take out zeroed terms
call psb_spcnv(af,info,afmt='coo')
call acsrf%mv_to_coo(acoof,info)
k = 0
do j=1,psb_sp_get_nnzeros(af)
if ((af%aspk(j) /= dzero) .or. (af%ia1(j)==af%ia2(j))) then
do j=1,acoof%get_nzeros()
if ((acoof%val(j) /= dzero) .or. (acoof%ia(j)==acoof%ja(j))) then
k = k + 1
af%aspk(k) = af%aspk(j)
af%ia1(k) = af%ia1(j)
af%ia2(k) = af%ia2(j)
acoof%val(k) = acoof%val(j)
acoof%ia(k) = acoof%ia(j)
acoof%ja(k) = acoof%ja(j)
end if
end do
!!$ write(debug_unit,*) me,' ',trim(name),' Non zeros from filtered matrix:',k,af%m,af%k
call psb_sp_setifld(k,psb_nnz_,af,info)
call psb_spcnv(af,info,afmt='csr')
call acoof%set_nzeros(k)
call acsrf%mv_from_coo(acoof,info)
end if
@ -301,9 +277,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
end do
if (filter_mat) call psb_sp_scal(adiag,af,info)
call psb_sp_scal(adiag,am3,info)
if (filter_mat) call acsrf%scal(adiag,info)
if (info == 0) call acsr3%scal(adiag,info)
if (info /= 0) goto 9999
@ -316,30 +291,25 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! This only works with CSR
!
if (psb_toupper(am3%fida)=='CSR') then
anorm = dzero
dg = done
do i=1,am3%m
tmp = dzero
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) <= am3%m) then
tmp = tmp + dabs(am3%aspk(j))
endif
if (am3%ia1(j) == i ) then
dg = dabs(am3%aspk(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
info = 4001
call psb_errpush(info,name,a_err='this section only CSR')
goto 9999
endif
anorm = dzero
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = dzero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + dabs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = dabs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
anorm = psb_spnrmi(am3,desc_a,info)
anorm = acsr3%csnmi()
endif
if (info /= 0) then
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
@ -368,20 +338,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Build the smoothed prolongator using the filtered matrix
!
if (psb_toupper(af%fida)=='CSR') then
do i=1,af%m
do j=af%ia2(i),af%ia2(i+1)-1
if (af%ia1(j) == i) then
af%aspk(j) = done - omega*af%aspk(j)
else
af%aspk(j) = - omega*af%aspk(j)
end if
end do
do i=1,acsrf%get_nrows()
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) then
acsrf%val(j) = done - omega*acsrf%val(j)
else
acsrf%val(j) = - omega*acsrf%val(j)
end if
end do
else
call psb_errpush(4001,name,a_err='Invalid AF storage format')
goto 9999
end if
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -389,17 +354,17 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Symbmm90 does the allocation for its result.
!
! am1 = (I-w*D*Af)Ptilde
! acsrm1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(af,am4,am1,info)
call psb_symbmm(acsrf,acsr4,acsr1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(af,am4,am1)
call psb_numbmm(acsrf,acsr4,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -408,20 +373,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Build the smoothed prolongator using the original matrix
!
if (psb_toupper(am3%fida)=='CSR') then
do i=1,am3%m
do j=am3%ia2(i),am3%ia2(i+1)-1
if (am3%ia1(j) == i) then
am3%aspk(j) = done - omega*am3%aspk(j)
else
am3%aspk(j) = - omega*am3%aspk(j)
end if
end do
do i=1,acsr3%get_nrows()
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) == i) then
acsr3%val(j) = done - omega*acsr3%val(j)
else
acsr3%val(j) = - omega*acsr3%val(j)
end if
end do
else
call psb_errpush(4001,name,a_err='Invalid AM3 storage format')
goto 9999
end if
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -429,30 +389,26 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Symbmm90 does the allocation for its result.
!
! am1 = (I-w*D*A)Ptilde
! acsrm1 = (I-w*D*A)Ptilde
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(am3,am4,am1,info)
call psb_symbmm(acsr3,acsr4,acsr1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(am3,am4,am1)
call psb_numbmm(acsr3,acsr4,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
call acsr4%free()
call psb_sp_free(am4,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_free')
goto 9999
end if
call am1%mv_from(acsr1)
if (ml_global_nmb) then
!
! Now we have to gather the halo of am1, and add it to itself
@ -461,7 +417,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_sphalo(am1,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if (info == 0) call psb_rwextd(ncol,am1,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
if (info == 0) call am4%free()
else
call psb_rwextd(ncol,am1,info)
endif
@ -482,29 +438,32 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
& 'Done NUMBMM 2'
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
call psb_transp(am1,am2,fmt='COO')
nzl = am2%infoa(psb_nnz_)
call am2%transp(am1)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < am2%ia1(k)) .and.(am2%ia1(k) <= naggrp1)) then
if ((naggrm1 < acoo2%ia(k)) .and.(acoo2%ia(k) <= naggrp1)) then
i = i+1
am2%aspk(i) = am2%aspk(k)
am2%ia1(i) = am2%ia1(k)
am2%ia2(i) = am2%ia2(k)
acoo2%val(i) = acoo2%val(k)
acoo2%ia(i) = acoo2%ia(k)
acoo2%ja(i) = acoo2%ja(k)
end if
end do
am2%infoa(psb_nnz_) = i
call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
call acoo2%set_nzeros(i)
call acoo2%trim()
call am2%mv_from(acoo2)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /=0) then
call psb_errpush(4010,name,a_err='spcnv am2')
goto 9999
end if
else
call psb_transp(am1,am2)
call am2%transp(am1)
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -515,7 +474,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if (info == 0) call psb_rwextd(ncol,am3,info,b=am4)
if (info == 0) call psb_sp_free(am4,info)
if (info == 0) call am4%free()
else if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
endif
@ -530,8 +489,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
& 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info)
if (info == 0) call psb_numbmm(am2,am3,b)
if (info == 0) call psb_sp_free(am3,info)
if (info == 0) call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (info == 0) call am3%free()
if (info == 0) call b%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4001,name,a_err='Build b = am2 x am3')
goto 9999
@ -547,14 +506,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,p%ac,info)
nzac = p%ac%infoa(psb_nnz_)
nzl = p%ac%infoa(psb_nnz_)
nzac = b%get_nzeros()
nzl = nzac
call b%mv_to(bcoo)
if (info == 0) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == 0) call psb_cdins(nzl,p%ac%ia1,p%ac%ia2,p%desc_ac,info)
if (info == 0) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == 0) call psb_cdasb(p%desc_ac,info)
if (info == 0) call psb_glob_to_loc(p%ac%ia1(1:nzl),p%desc_ac,info,iact='I')
if (info == 0) call psb_glob_to_loc(p%ac%ia2(1:nzl),p%desc_ac,info,iact='I')
if (info == 0) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == 0) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= 0) then
call psb_errpush(4001,name,a_err='Creating p%desc_ac and converting ac')
goto 9999
@ -562,14 +522,12 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
call p%ac%set_nrows(psb_cd_get_local_rows(p%desc_ac))
call p%ac%set_ncols(psb_cd_get_local_cols(p%desc_ac))
call p%ac%set_asb()
p%ac%m=psb_cd_get_local_rows(p%desc_ac)
p%ac%k=psb_cd_get_local_cols(p%desc_ac)
p%ac%fida='COO'
p%ac%descra='GUN'
call psb_sp_free(b,info)
if (info == 0) deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
@ -577,26 +535,31 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
if (np>1) then
nzl = psb_sp_get_nnzeros(am1)
call psb_glob_to_loc(am1%ia1(1:nzl),p%desc_ac,info,'I')
call am1%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_glob_to_loc')
goto 9999
end if
call am1%mv_from(acsr1)
endif
am1%k=psb_cd_get_local_cols(p%desc_ac)
call am1%set_ncols(psb_cd_get_local_cols(p%desc_ac))
if (np>1) then
call psb_spcnv(am2,info,afmt='coo',dupl=psb_dupl_add_)
nzl = am2%infoa(psb_nnz_)
if (info == 0) call psb_glob_to_loc(am2%ia1(1:nzl),p%desc_ac,info,'I')
if (info == 0) call psb_spcnv(am2,info,afmt='csr',dupl=psb_dupl_add_)
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
if (info == 0) call psb_glob_to_loc(acoo2%ia(1:nzl),p%desc_ac,info,'I')
call acoo2%set_dupl(psb_dupl_add_)
if (info == 0) call am2%mv_from(acoo2)
if (info == 0) call am2%cscnv(info,type='csr')
if(info /= 0) then
call psb_errpush(4001,name,a_err='Converting am2 to local')
goto 9999
end if
end if
am2%m=psb_cd_get_local_cols(p%desc_ac)
call am2%set_nrows(psb_cd_get_local_cols(p%desc_ac))
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -606,39 +569,43 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == 0) call psb_sp_all(ntaggr,ntaggr,p%ac,nzac,info)
if (info /= 0) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,p%ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,p%ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,p%ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
if (.false.) then
nzbr(:) = 0
nzbr(me+1) = b%get_nzeros()
call b%mv_to(bcoo)
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == 0) call cootmp%allocate(ntaggr,ntaggr,nzac)
if (info /= 0) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,&
& cootmp%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == 0) call mpi_allgatherv(bcoo%ia,ndx,mpi_integer,&
& cootmp%ia,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info == 0) call mpi_allgatherv(bcoo%ja,ndx,mpi_integer,&
& cootmp%ja,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
call bcoo%free()
p%ac%m = ntaggr
p%ac%k = ntaggr
p%ac%infoa(psb_nnz_) = nzac
p%ac%fida='COO'
p%ac%descra='GUN'
call psb_spcnv(p%ac,info,afmt='coo',dupl=psb_dupl_add_)
call cootmp%set_nzeros(nzac)
call cootmp%set_dupl(psb_dupl_add_)
call p%ac%mv_from(cootmp)
else
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
endif
if(info /= 0) goto 9999
call psb_sp_free(b,info)
if(info /= 0) goto 9999
deallocate(nzbr,idisp,stat=info)
@ -660,10 +627,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_)
call psb_sp_clone(b,p%ac,info)
call psb_move_alloc(b,p%ac,info)
if (info == 0) call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == 0) call psb_cdasb(p%desc_ac,info)
if (info == 0) call psb_sp_free(b,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Build desc_ac, ac')
goto 9999
@ -678,47 +644,51 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(4010,name,a_err='psb_cdall')
goto 9999
end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= 0) goto 9999
!!$
!!$ nzbr(:) = 0
!!$ nzbr(me+1) = b%get_nzeros()
!!$ call psb_sum(ictxt,nzbr(1:np))
!!$ nzac = sum(nzbr)
!!$
!!$ call b%mv_to(bcoo)
!!$ call psb_sum(ictxt,nzbr(1:np))
!!$ nzac = sum(nzbr)
!!$ if (info == 0) call cootmp%allocate(ntaggr,ntaggr,nzac)
!!$ if (info /= 0) goto 9999
!!$
!!$ do ip=1,np
!!$ idisp(ip) = sum(nzbr(1:ip-1))
!!$ enddo
!!$ ndx = nzbr(me+1)
!!$
!!$ call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,&
!!$ & cootmp%val,nzbr,idisp,&
!!$ & mpi_double_precision,icomm,info)
!!$ if (info == 0) call mpi_allgatherv(bcoo%ia,ndx,mpi_integer,&
!!$ & cootmp%ia,nzbr,idisp,&
!!$ & mpi_integer,icomm,info)
!!$
!!$ if (info == 0) call mpi_allgatherv(bcoo%ja,ndx,mpi_integer,&
!!$ & cootmp%ja,nzbr,idisp,&
!!$ & mpi_integer,icomm,info)
!!$
!!$ if (info /= 0) then
!!$ call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
!!$ goto 9999
!!$ end if
!!$ call bcoo%free()
!!$ call cootmp%set_nzeros(nzac)
!!$ call cootmp%set_dupl(psb_dupl_add_)
!!$ call p%ac%mv_from(cootmp)
!!$ if(info /= 0) goto 9999
nzbr(:) = 0
nzbr(me+1) = b%infoa(psb_nnz_)
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call psb_sp_all(ntaggr,ntaggr,p%ac,nzac,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_all')
goto 9999
end if
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(b%aspk,ndx,mpi_double_precision,p%ac%aspk,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == 0) call mpi_allgatherv(b%ia1,ndx,mpi_integer,p%ac%ia1,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info == 0) call mpi_allgatherv(b%ia2,ndx,mpi_integer,p%ac%ia2,nzbr,idisp,&
& mpi_integer,icomm,info)
deallocate(nzbr,idisp,stat=info)
if (info /= 0) then
call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
goto 9999
end if
p%ac%m = ntaggr
p%ac%k = ntaggr
p%ac%infoa(psb_nnz_) = nzac
p%ac%fida='COO'
p%ac%descra='GUN'
call psb_spcnv(p%ac,info,afmt='coo',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
end if
call psb_sp_free(b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
info = 4000
call psb_errpush(info,name)
goto 9999
end if
@ -742,7 +712,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end select
call psb_spcnv(p%ac,info,afmt='csr',dupl=psb_dupl_add_)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spcnv')
goto 9999
@ -755,8 +725,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr)
if (info == 0) call psb_sp_free(am1,info)
if (info == 0) call psb_sp_free(am2,info)
if (info == 0) call am1%free()
if (info == 0) call am2%free()
if(info /= 0) then
call psb_errpush(4010,name,a_err='sp_Free')
goto 9999

12
mlprec/mld_das_aply.f90
Unescape Escape View File

@ -203,7 +203,8 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! local extended matrix, stored into the permutation vector prec%perm
!
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%perm,tx,info)
!!$ call psb_gelp('n',prec%perm,tx,info)
info = 1
if(info /=0) then
info=4010
ch_err='psb_gelp'
@ -227,7 +228,8 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Apply to ty the inverse permutation of prec%perm
!
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%invperm,ty,info)
!!$ call psb_gelp('n',prec%invperm,ty,info)
info = 1
if(info /= 0) then
info=4010
ch_err='psb_gelp'
@ -313,7 +315,8 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! local extended matrix, stored into the permutation vector prec%perm
!
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%perm,tx,info)
!!$ call psb_gelp('n',prec%perm,tx,info)
info = 1
if(info /=0) then
info=4010
ch_err='psb_gelp'
@ -337,7 +340,8 @@ subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Apply to ty the inverse permutation of prec%perm
!
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%invperm,ty,info)
!!$ call psb_gelp('n',prec%invperm,ty,info)
info = 1
if(info /= 0) then
info=4010
ch_err='psb_gelp'

45
mlprec/mld_das_bld.f90
Unescape Escape View File

@ -74,18 +74,18 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
Implicit None
! Arguments
type(psb_dspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
type(psb_d_sparse_mat), intent(in), target :: a
Type(psb_desc_type), Intent(in) :: desc_a
type(mld_dbaseprec_type), intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
character, intent(in) :: upd
integer, intent(out) :: info
! Local variables
integer :: ptype,novr
integer :: icomm
Integer :: np,me,nnzero,ictxt, int_err(5),&
& tot_recv, n_row,n_col,nhalo, err_act, data_
type(psb_dspmat_type) :: blck
type(psb_d_sparse_mat) :: blck
integer :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -108,7 +108,7 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
nnzero = psb_sp_get_nnzeros(a)
nnzero = a%get_nzeros()
nhalo = n_col-n_row
ptype = p%iprcparm(mld_smoother_type_)
novr = p%iprcparm(mld_sub_ovr_)
@ -138,17 +138,8 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
& write(debug_unit,*) me,' ',trim(name),&
& 'Early return: P>=3 N_OVR=0'
endif
call psb_sp_all(0,0,blck,1,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
blck%fida = 'COO'
blck%infoa(psb_nnz_) = 0
call mld_fact_bld(a,p,upd,info,blck=blck)
call mld_fact_bld(a,p,upd,info)
if (info /= 0) then
info=4010
@ -192,15 +183,7 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
& write(debug_unit,*) me,' ',trim(name),&
& 'Early return: P>=3 N_OVR=0'
endif
call psb_sp_all(0,0,blck,1,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
blck%fida = 'COO'
blck%infoa(psb_nnz_) = 0
call blck%csall(0,0,info,1)
else
@ -216,8 +199,8 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
call psb_cdbldext(a,desc_a,novr,p%desc_data,info,extype=psb_ovt_asov_)
if(debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' From cdbldext _:',p%desc_data%matrix_data(psb_n_row_),&
& p%desc_data%matrix_data(psb_n_col_)
& ' From cdbldext _:',psb_cd_get_local_rows(p%desc_data),&
& psb_cd_get_local_cols(p%desc_data)
if (info /= 0) then
info=4010
@ -229,7 +212,7 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Before sphalo ',blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_)
& 'Before sphalo '
!
! Retrieve the remote sparse matrix rows required for the AS extended
@ -244,10 +227,10 @@ subroutine mld_das_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
if (debug_level >=psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'After psb_sphalo ',&
& blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_)
& blck%get_nrows(), blck%get_nzeros()
End if

4
mlprec/mld_dbaseprec_bld.f90
Unescape Escape View File

@ -76,9 +76,9 @@ subroutine mld_dbaseprec_bld(a,desc_a,p,info,upd)
Implicit None
! Arguments
type(psb_dspmat_type), target :: a
type(psb_d_sparse_mat), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dbaseprec_type),intent(inout) :: p
type(mld_dbaseprec_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd

14
mlprec/mld_dcoarse_bld.f90
Unescape Escape View File

@ -73,17 +73,15 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_donelev_type), intent(inout),target :: p
integer, intent(out) :: info
integer, intent(out) :: info
! Local variables
type(psb_desc_type) :: desc_ac
type(psb_dspmat_type) :: ac
character(len=20) :: name
integer :: ictxt, np, me, err_act
integer, allocatable :: ilaggr(:), nlaggr(:)
character(len=20) :: name
integer :: ictxt, np, me, err_act
integer, allocatable :: ilaggr(:), nlaggr(:)
name='mld_dcoarse_bld'
if (psb_get_errstatus().ne.0) return

40
mlprec/mld_ddiag_bld.f90
Unescape Escape View File

@ -65,10 +65,10 @@ subroutine mld_ddiag_bld(a,desc_a,p,info)
Implicit None
! Arguments
type(psb_dspmat_type),intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(mld_dbaseprec_type),intent(inout) :: p
integer, intent(out) :: info
type(psb_d_sparse_mat),intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(mld_dbaseprec_type),intent(inout) :: p
integer, intent(out) :: info
! Local variables
Integer :: err_act,ictxt, me, np, n_row, n_col,i
@ -99,7 +99,7 @@ subroutine mld_ddiag_bld(a,desc_a,p,info)
!
! Retrieve the diagonal entries of the matrix A
!
call psb_sp_getdiag(a,p%d,info)
call a%get_diag(p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_getdiag'
@ -128,21 +128,21 @@ subroutine mld_ddiag_bld(a,desc_a,p,info)
endif
end do
if (a%pl(1) /= 0) then
!
! Apply the same row permutation as in the sparse matrix A
!
call psb_gelp('n',a%pl,p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_gelp'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
endif
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),'Done'
!!$ if (a%pl(1) /= 0) then
!!$ !
!!$ ! Apply the same row permutation as in the sparse matrix A
!!$ !
!!$ call psb_gelp('n',a%pl,p%d,info)
!!$ if(info /= 0) then
!!$ info=4010
!!$ ch_err='psb_gelp'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ end if
!!$ endif
!!$
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),'Done'
call psb_erractionrestore(err_act)
return

106
mlprec/mld_dfact_bld.f90
Unescape Escape View File

@ -116,20 +116,20 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
implicit none
! Arguments
type(psb_dspmat_type), intent(in), target :: a
type(psb_d_sparse_mat), intent(in), target :: a
type(mld_dbaseprec_type), intent(inout) :: p
integer, intent(out) :: info
character, intent(in) :: upd
type(psb_dspmat_type), intent(in), target, optional :: blck
type(psb_d_sparse_mat), intent(in), target, optional :: blck
! Local Variables
type(psb_dspmat_type), pointer :: blck_
type(psb_dspmat_type) :: atmp
integer :: ictxt,np,me,err_act
integer :: debug_level, debug_unit
integer :: k, m, int_err(5), n_row, nrow_a, n_col
character :: trans, unitd
character(len=20) :: name, ch_err
type(psb_d_sparse_mat), pointer :: blck_
type(psb_d_sparse_mat) :: atmp
integer :: ictxt,np,me,err_act
integer :: debug_level, debug_unit
integer :: k, m, int_err(5), n_row, nrow_a, n_col
character :: trans, unitd
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
@ -140,7 +140,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
ictxt = psb_cd_get_context(p%desc_data)
call psb_info(ictxt, me, np)
m = a%m
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
@ -155,17 +155,15 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
blck_ => blck
else
allocate(blck_,stat=info)
if (info ==0) call psb_sp_all(0,0,blck_,1,info)
if (info ==0) call blck_%csall(0,0,info,1)
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
blck_%fida = 'COO'
blck_%infoa(psb_nnz_) = 0
end if
call psb_nullify_sp(atmp)
!
! Treat separately the case the local matrix has to be reordered
@ -194,16 +192,16 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
! matrix. The clipped matrix is then stored in CSR format.
!
if (p%iprcparm(mld_smoother_sweeps_) > 1) then
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
call atmp%csclip(p%av(mld_ap_nd_),info,&
& jmin=atmp%get_nrows()+1,rscale=.false.,cscale=.false.)
if (info == 0) call p%av(mld_ap_nd_)%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
k = p%av(mld_ap_nd_)%get_nzeros()
call psb_sum(ictxt,k)
if (k == 0) then
@ -215,9 +213,6 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
p%iprcparm(mld_smoother_sweeps_) = 1
end if
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' Factoring rows ',&
& atmp%m,a%m,blck_%m,atmp%ia2(atmp%m+1)-1
!
! Compute a factorization of the diagonal block of the local matrix,
@ -229,7 +224,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
! ILU(k)/MILU(k)/ILU(k,t) factorization.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_ilu_bld(atmp,p,upd,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='mld_ilu_bld')
@ -240,7 +235,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
! LU factorization through the SuperLU package.
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
@ -253,7 +248,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
! when the matrix is distributed among the processes.
! NOTE: Should have NO overlap here!!!!
!
call psb_spcnv(a,atmp,info,afmt='csr')
call a%cscnv(atmp,info,type='csr')
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_sludist_bld')
@ -264,7 +259,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
! LU factorization through the UMFPACK package.
!
call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_umf_bld')
@ -285,12 +280,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
goto 9999
end select
call psb_sp_free(atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call atmp%free()
!
! No reordering of the local matrix is required
@ -305,22 +295,22 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
if (p%iprcparm(mld_smoother_sweeps_) > 1) then
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
nrow_a = a%m
nrow_a = a%get_nrows()
! The following is known to work
! given that the output from CLIP is in COO.
call psb_sp_clip(a,p%av(mld_ap_nd_),info,&
call a%csclip(p%av(mld_ap_nd_),info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
if (info == 0) call psb_sp_clip(blck_,atmp,info,&
if (info == 0) call blck_%csclip(atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
if (info == 0) call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,&
& afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call p%av(mld_ap_nd_)%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
k = p%av(mld_ap_nd_)%get_nzeros()
call psb_sum(ictxt,k)
if (k == 0) then
@ -331,11 +321,8 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
p%iprcparm(mld_smoother_sweeps_) = 1
end if
call psb_sp_free(atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call atmp%free()
end if
!
! Compute a factorization of the diagonal block of the local matrix,
@ -362,24 +349,21 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_spcnv(a,atmp,info,afmt='coo')
call a%cscnv(atmp,info,type='coo')
if (info == 0) call psb_rwextd(n_row,atmp,info,b=blck_)
!
! Compute the LU factorization.
!
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call atmp%free()
case(mld_sludist_)
!
@ -387,25 +371,21 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
! when the matrix is distributed among the processes.
! NOTE: Should have NO overlap here!!!!
!
call psb_spcnv(a,atmp,info,afmt='csr')
call a%cscnv(atmp,info,type='csr')
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='mld_sludist_bld')
goto 9999
end if
call psb_sp_free(atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call atmp%free()
case(mld_umf_)
!
! LU factorization through the UMFPACK package.
!
call psb_spcnv(a,atmp,info,afmt='coo')
call a%cscnv(atmp,info,type='coo')
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv')
goto 9999
@ -418,7 +398,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
!
! Compute the LU factorization.
!
if (info == 0) call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
if (info == 0) call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -428,11 +408,7 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
goto 9999
end if
call psb_sp_free(atmp,info)
if (info/=0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999
end if
call atmp%free()
case(mld_f_none_)
!
@ -455,8 +431,8 @@ subroutine mld_dfact_bld(a,p,upd,info,blck)
end select
if (.not.present(blck)) then
call psb_sp_free(blck_,info)
if (info == 0) deallocate(blck_)
call blck_%free()
if (info == 0) deallocate(blck_, stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_sp_free')
goto 9999

532
mlprec/mld_dilu0_fact.f90
Unescape Escape View File

@ -99,7 +99,7 @@
! greater than 0. If the overlap is 0 or the matrix has been reordered
! (see mld_fact_bld), then blck is empty.
!
subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck)
subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck, upd)
use psb_base_mod
use mld_inner_mod, mld_protect_name => mld_dilu0_fact
@ -107,17 +107,20 @@ subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck)
implicit none
! Arguments
integer, intent(in) :: ialg
type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u
real(psb_dpk_), intent(inout) :: d(:)
integer, intent(out) :: info
type(psb_dspmat_type),intent(in), optional, target :: blck
integer, intent(in) :: ialg
type(psb_d_sparse_mat),intent(in) :: a
type(psb_d_sparse_mat),intent(inout) :: l,u
real(psb_dpk_), intent(inout) :: d(:)
integer, intent(out) :: info
type(psb_d_sparse_mat),intent(in), optional, target :: blck
character, intent(in), optional :: upd
! Local variables
integer :: l1, l2,m,err_act
type(psb_dspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
integer :: l1, l2, m, err_act
type(psb_d_sparse_mat), pointer :: blck_
type(psb_d_csr_sparse_mat) :: ll, uu
character :: upd_
character(len=20) :: name, ch_err
name='mld_dilu0_fact'
info = 0
@ -130,28 +133,36 @@ subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck)
blck_ => blck
else
allocate(blck_,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
call psb_nullify_sp(blck_) ! Probably pointless.
call psb_sp_all(0,0,blck_,1,info)
if(info.ne.0) then
if (info == 0) call blck_%csall(0,0,info,1)
if (info /= 0) then
info=4010
ch_err='psb_sp_all'
ch_err='csall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
blck_%m=0
endif
if (present(upd)) then
upd_ = psb_toupper(upd)
else
upd_ = 'F'
end if
m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
& (m > size(d)) ) then
write(0,*) 'Wrong allocation status for L,D,U? ',&
& l%get_nrows(),size(d),u%get_nrows()
info = -1
return
end if
call l%mv_to(ll)
call u%mv_to(uu)
!
! Compute the ILU(0) or the MILU(0) factorization, depending on ialg
!
call mld_dilu0_factint(ialg,m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
call mld_dilu0_factint(ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,upd_,info)
if(info.ne.0) then
info=4010
ch_err='mld_dilu0_factint'
@ -162,24 +173,22 @@ subroutine mld_dilu0_fact(ialg,a,l,u,d,info,blck)
!
! Store information on the L and U sparse matrices
!
l%infoa(1) = l1
l%fida = 'CSR'
l%descra = 'TLU'
u%infoa(1) = l2
u%fida = 'CSR'
u%descra = 'TUU'
l%m = m
l%k = m
u%m = m
u%k = m
call l%mv_from(ll)
call l%set_triangle()
call l%set_unit()
call l%set_lower()
call u%mv_from(uu)
call u%set_triangle()
call u%set_unit()
call u%set_upper()
!
! Nullify pointer / deallocate memory
!
if (present(blck)) then
blck_ => null()
else
call psb_sp_free(blck_,info)
call blck_%free()
if(info.ne.0) then
info=4010
ch_err='psb_sp_free'
@ -220,7 +229,7 @@ contains
! solve stage associated to the ILU(0)/MILU(0) factorization).
!
! The routine copies and factors "on the fly" from the sparse matrix structures a
! and b into the arrays laspk, uaspk, d (L, U without its diagonal, diagonal of U).
! and b into the arrays lval, uval, d (L, U without its diagonal, diagonal of U).
!
!
! Arguments:
@ -252,49 +261,50 @@ contains
! d - real(psb_dpk_), dimension(:), output.
! The inverse of the diagonal entries of the U factor in the
! incomplete factorization.
! laspk - real(psb_dpk_), dimension(:), input/output.
! lval - real(psb_dpk_), dimension(:), input/output.
! The entries of U are stored according to the CSR format.
! The L factor in the incomplete factorization.
! lia1 - integer, dimension(:), input/output.
! lja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the L factor,
! according to the CSR storage format.
! lia2 - integer, dimension(:), input/output.
! lirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the L factor in laspk, according to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! of the L factor in lval, according to the CSR storage format.
! uval - real(psb_dpk_), dimension(:), input/output.
! The U factor in the incomplete factorization.
! The entries of U are stored according to the CSR format.
! uia1 - integer, dimension(:), input/output.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the U factor,
! according to the CSR storage format.
! uia2 - integer, dimension(:), input/output.
! uirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the U factor in uaspk, according to the CSR storage format.
! of the U factor in uval, according to the CSR storage format.
! l1 - integer, output.
! The number of nonzero entries in laspk.
! The number of nonzero entries in lval.
! l2 - integer, output.
! The number of nonzero entries in uaspk.
! The number of nonzero entries in uval.
! info - integer, output.
! Error code.
!
subroutine mld_dilu0_factint(ialg,m,ma,a,mb,b,&
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
subroutine mld_dilu0_factint(ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info)
implicit none
! Arguments
integer, intent(in) :: ialg
type(psb_dspmat_type),intent(in) :: a,b
integer,intent(inout) :: m,l1,l2,info
integer, intent(in) :: ma,mb
integer, dimension(:), intent(inout) :: lia1,lia2,uia1,uia2
real(psb_dpk_), dimension(:),intent(inout) :: laspk,uaspk,d
integer, intent(in) :: ialg
type(psb_d_sparse_mat),intent(in) :: a,b
integer,intent(inout) :: l1,l2,info
integer, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_dpk_), intent(inout) :: lval(:),uval(:),d(:)
character, intent(in) :: upd
! Local variables
integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act
integer :: i,j,k,l,low1,low2,kk,jj,ll, ktrw,err_act, m
integer :: ma,mb
real(psb_dpk_) :: dia,temp
integer, parameter :: nrb=16
type(psb_dspmat_type) :: trw
type(psb_d_coo_sparse_mat) :: trw
integer :: int_err(5)
character(len=20) :: name, ch_err
@ -302,6 +312,8 @@ contains
if(psb_get_errstatus().ne.0) return
info=0
call psb_erractionsave(err_act)
ma = a%get_nrows()
mb = b%get_nrows()
select case(ialg)
case(mld_ilu_n_,mld_milu_n_)
@ -312,154 +324,152 @@ contains
goto 9999
end select
call psb_nullify_sp(trw)
trw%m=0
trw%k=0
call psb_sp_all(trw,1,info)
if(info.ne.0) then
call trw%allocate(0,0,1)
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
lia2(1) = 1
uia2(1) = 1
l1 = 0
l2 = 0
m = ma+mb
!
! Cycle over the matrix rows
!
do i = 1, m
if (psb_toupper(upd) == 'F' ) then
lirp(1) = 1
uirp(1) = 1
l1 = 0
l2 = 0
!
! Cycle over the matrix rows
!
do i = 1, m
d(i) = dzero
d(i) = dzero
if (i <= ma) then
!
! Copy the i-th local row of the matrix, stored in a,
! into laspk/d(i)/uaspk
!
call ilu_copyin(i,ma,a,i,1,m,l1,lia1,laspk,&
& d(i),l2,uia1,uaspk,ktrw,trw)
else
!
! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into laspk/d(i)/uaspk
!
call ilu_copyin(i-ma,mb,b,i,1,m,l1,lia1,laspk,&
& d(i),l2,uia1,uaspk,ktrw,trw)
endif
if (i <= ma) then
!
! Copy the i-th local row of the matrix, stored in a,
! into lval/d(i)/uval
!
call ilu_copyin(i,ma,a,i,1,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd)
else
!
! Copy the i-th local row of the matrix, stored in b
! (as (i-ma)-th row), into lval/d(i)/uval
!
call ilu_copyin(i-ma,mb,b,i,1,m,l1,lja,lval,&
& d(i),l2,uja,uval,ktrw,trw,upd)
endif
lia2(i+1) = l1 + 1
uia2(i+1) = l2 + 1
lirp(i+1) = l1 + 1
uirp(i+1) = l2 + 1
dia = d(i)
do kk = lia2(i), lia2(i+1) - 1
!
! Compute entry l(i,k) (lower factor L) of the incomplete
! factorization
!
temp = laspk(kk)
k = lia1(kk)
laspk(kk) = temp*d(k)
!
! Update the rest of row i (lower and upper factors L and U)
! using l(i,k)
!
low1 = kk + 1
low2 = uia2(i)
!
updateloop: do jj = uia2(k), uia2(k+1) - 1
dia = d(i)
do kk = lirp(i), lirp(i+1) - 1
!
j = uia1(jj)
! Compute entry l(i,k) (lower factor L) of the incomplete
! factorization
!
if (j < i) then
!
! search l(i,*) (i-th row of L) for a matching index j
!
do ll = low1, lia2(i+1) - 1
l = lia1(ll)
if (l > j) then
low1 = ll
exit
else if (l == j) then
laspk(ll) = laspk(ll) - temp*uaspk(jj)
low1 = ll + 1
cycle updateloop
end if
enddo
else if (j == i) then
temp = lval(kk)
k = lja(kk)
lval(kk) = temp*d(k)
!
! Update the rest of row i (lower and upper factors L and U)
! using l(i,k)
!
low1 = kk + 1
low2 = uirp(i)
!
updateloop: do jj = uirp(k), uirp(k+1) - 1
!
! j=i: update the diagonal
j = uja(jj)
!
dia = dia - temp*uaspk(jj)
cycle updateloop
if (j < i) then
!
! search l(i,*) (i-th row of L) for a matching index j
!
do ll = low1, lirp(i+1) - 1
l = lja(ll)
if (l > j) then
low1 = ll
exit
else if (l == j) then
lval(ll) = lval(ll) - temp*uval(jj)
low1 = ll + 1
cycle updateloop
end if
enddo
else if (j == i) then
!
! j=i: update the diagonal
!
dia = dia - temp*uval(jj)
cycle updateloop
!
else if (j > i) then
!
! search u(i,*) (i-th row of U) for a matching index j
!
do ll = low2, uirp(i+1) - 1
l = uja(ll)
if (l > j) then
low2 = ll
exit
else if (l == j) then
uval(ll) = uval(ll) - temp*uval(jj)
low2 = ll + 1
cycle updateloop
end if
enddo
end if
!
else if (j > i) then
!
! search u(i,*) (i-th row of U) for a matching index j
! If we get here we missed the cycle updateloop, which means
! that this entry does not match; thus we accumulate on the
! diagonal for MILU(0).
!
do ll = low2, uia2(i+1) - 1
l = uia1(ll)
if (l > j) then
low2 = ll
exit
else if (l == j) then
uaspk(ll) = uaspk(ll) - temp*uaspk(jj)
low2 = ll + 1
cycle updateloop
end if
enddo
end if
if (ialg == mld_milu_n_) then
dia = dia - temp*uval(jj)
end if
enddo updateloop
enddo
!
! Check the pivot size
!
if (abs(dia) < epstol) then
!
! Too small pivot: unstable factorization
!
! If we get here we missed the cycle updateloop, which means
! that this entry does not match; thus we accumulate on the
! diagonal for MILU(0).
info = 2
int_err(1) = i
write(ch_err,'(g20.10)') abs(dia)
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
goto 9999
else
!
if (ialg == mld_milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloop
enddo
!
! Check the pivot size
!
if (abs(dia) < epstol) then
!
! Too small pivot: unstable factorization
!
info = 2
int_err(1) = i
write(ch_err,'(g20.10)') abs(dia)
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
goto 9999
else
! Compute 1/pivot
!
dia = done/dia
end if
d(i) = dia
!
! Compute 1/pivot
! Scale row i of upper triangle
!
dia = done/dia
end if
d(i) = dia
!
! Scale row i of upper triangle
!
do kk = uia2(i), uia2(i+1) - 1
uaspk(kk) = uaspk(kk)*dia
do kk = uirp(i), uirp(i+1) - 1
uval(kk) = uval(kk)*dia
enddo
enddo
enddo
call psb_sp_free(trw,info)
if(info.ne.0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
else
write(0,*) 'Update not implemented '
info = 31
call psb_errpush(info,name,i_err=(/13,0,0,0,0/),a_err=upd)
goto 9999
end if
call trw%free()
call psb_erractionrestore(err_act)
return
@ -478,15 +488,15 @@ contains
! Note: internal subroutine of mld_dilu0_fact
!
! This routine copies a row of a sparse matrix A, stored in the psb_dspmat_type
! data structure a, into the arrays laspk and uaspk and into the scalar variable
! data structure a, into the arrays lval and uval and into the scalar variable
! dia, corresponding to the lower and upper triangles of A and to the diagonal
! entry of the row, respectively. The entries in laspk and uaspk are stored
! entry of the row, respectively. The entries in lval and uval are stored
! according to the CSR format; the corresponding column indices are stored in
! the arrays lia1 and uia1.
! the arrays lja and uja.
!
! If the sparse matrix is in CSR format, a 'straight' copy is performed;
! otherwise psb_sp_getblk is used to extract a block of rows, which is then
! copied into laspk, dia, uaspk row by row, through successive calls to
! copied into lval, dia, uval row by row, through successive calls to
! ilu_copyin.
!
! The routine is used by mld_dilu0_factint in the computation of the ILU(0)/MILU(0)
@ -514,23 +524,23 @@ contains
! The output matrix will contain a clipped copy taken from
! a(1:m,jmin:jmax).
! l1 - integer, input/output.
! Pointer to the last occupied entry of laspk.
! lia1 - integer, dimension(:), input/output.
! Pointer to the last occupied entry of lval.
! lja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the lower triangle
! copied in laspk row by row (see mld_dilu0_factint), according
! copied in lval row by row (see mld_dilu0_factint), according
! to the CSR storage format.
! laspk - real(psb_dpk_), dimension(:), input/output.
! lval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! lower triangle are copied.
! dia - real(psb_dpk_), output.
! The diagonal entry of the copied row.
! l2 - integer, input/output.
! Pointer to the last occupied entry of uaspk.
! uia1 - integer, dimension(:), input/output.
! Pointer to the last occupied entry of uval.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the upper triangle
! copied in uaspk row by row (see mld_dilu0_factint), according
! copied in uval row by row (see mld_dilu0_factint), according
! to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! uval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! upper triangle are copied.
! ktrw - integer, input/output.
@ -544,92 +554,104 @@ contains
! until we empty the buffer. Thus we will make a call to psb_sp_getblk
! every nrb calls to copyin. If A is in CSR format it is unused.
!
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lia1,laspk,&
& dia,l2,uia1,uaspk,ktrw,trw)
subroutine ilu_copyin(i,m,a,jd,jmin,jmax,l1,lja,lval,&
& dia,l2,uja,uval,ktrw,trw,upd)
use psb_base_mod
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(inout) :: trw
type(psb_d_sparse_mat), intent(in) :: a
type(psb_d_coo_sparse_mat), intent(inout) :: trw
integer, intent(in) :: i,m,jd,jmin,jmax
integer, intent(inout) :: ktrw,l1,l2
integer, intent(inout) :: lia1(:), uia1(:)
real(psb_dpk_), intent(inout) :: laspk(:), uaspk(:), dia
integer, intent(inout) :: lja(:), uja(:)
real(psb_dpk_), intent(inout) :: lval(:), uval(:), dia
character, intent(in) :: upd
! Local variables
integer :: k,j,info,irb
integer, parameter :: nrb=16
integer :: k,j,info,irb, nz
integer, parameter :: nrb=40
character(len=20), parameter :: name='ilu_copyin'
character(len=20) :: ch_err
if (psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
if (psb_toupper(upd) == 'F') then
if (psb_toupper(a%fida)=='CSR') then
select type(aa => a%a)
type is (psb_d_csr_sparse_mat)
!
! Take a fast shortcut if the matrix is stored in CSR format
!
!
! Take a fast shortcut if the matrix is stored in CSR format
!
do j = a%ia2(i), a%ia2(i+1) - 1
k = a%ia1(j)
! write(0,*)'KKKKK',k
if ((k < jd).and.(k >= jmin)) then
l1 = l1 + 1
laspk(l1) = a%aspk(j)
lia1(l1) = k
else if (k == jd) then
dia = a%aspk(j)
else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1
uaspk(l2) = a%aspk(j)
uia1(l2) = k
end if
enddo
do j = aa%irp(i), aa%irp(i+1) - 1
k = aa%ja(j)
! write(0,*)'KKKKK',k
if ((k < jd).and.(k >= jmin)) then
l1 = l1 + 1
lval(l1) = aa%val(j)
lja(l1) = k
else if (k == jd) then
dia = aa%val(j)
else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1
uval(l2) = aa%val(j)
uja(l2) = k
end if
enddo
else
class default
!
! Otherwise use psb_sp_getblk, slower but able (in principle) of
! handling any format. In this case, a block of rows is extracted
! instead of a single row, for performance reasons, and these
! rows are copied one by one into laspk, dia, uaspk, through
! successive calls to ilu_copyin.
!
!
! Otherwise use psb_sp_getblk, slower but able (in principle) of
! handling any format. In this case, a block of rows is extracted
! instead of a single row, for performance reasons, and these
! rows are copied one by one into lval, dia, uval, through
! successive calls to ilu_copyin.
!
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if(info.ne.0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
ktrw=1
end if
do
if (ktrw > trw%infoa(psb_nnz_)) exit
if (trw%ia1(ktrw) > i) exit
k = trw%ia2(ktrw)
if ((k < jd).and.(k >= jmin)) then
l1 = l1 + 1
laspk(l1) = trw%aspk(ktrw)
lia1(l1) = k
else if (k == jd) then
dia = trw%aspk(ktrw)
else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1
uaspk(l2) = trw%aspk(ktrw)
uia1(l2) = k
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call aa%csget(i,i+irb-1,trw,info)
if(info /= 0) then
info=4010
ch_err='csget'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
ktrw=1
end if
ktrw = ktrw + 1
enddo
nz = trw%get_nzeros()
do
if (ktrw > nz) exit
if (trw%ia(ktrw) > i) exit
k = trw%ja(ktrw)
if ((k < jd).and.(k >= jmin)) then
l1 = l1 + 1
lval(l1) = trw%val(ktrw)
lja(l1) = k
else if (k == jd) then
dia = trw%val(ktrw)
else if ((k > jd).and.(k <= jmax)) then
l2 = l2 + 1
uval(l2) = trw%val(ktrw)
uja(l2) = k
end if
ktrw = ktrw + 1
enddo
end select
else
write(0,*) 'Update not implemented '
info = 31
call psb_errpush(info,name,i_err=(/13,0,0,0,0/),a_err=upd)
goto 9999
end if

251
mlprec/mld_dilu_bld.f90
Unescape Escape View File

@ -95,13 +95,13 @@ subroutine mld_dilu_bld(a,p,upd,info,blck)
use mld_inner_mod, mld_protect_name => mld_dilu_bld
implicit none
! Arguments
type(psb_dspmat_type), intent(in), target :: a
! Arguments
type(psb_d_sparse_mat), intent(in), target :: a
type(mld_dbaseprec_type), intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_dspmat_type), intent(in), optional :: blck
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_d_sparse_mat), intent(in), optional :: blck
! Local Variables
integer :: i, nztota, err_act, n_row, nrow_a
@ -123,143 +123,148 @@ subroutine mld_dilu_bld(a,p,upd,info,blck)
trans = 'N'
unitd = 'U'
!
! Check the memory available to hold the incomplete L and U factors
! and allocate it if needed
!
if (allocated(p%av)) then
if (size(p%av) < mld_bp_ilu_avsz_) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this. Just let it go.
! return
end if
enddo
deallocate(p%av,stat=info)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
n_row = psb_cd_get_local_rows(p%desc_data)
if (psb_toupper(upd) == 'F') then
!
! Check the memory available to hold the incomplete L and U factors
! and allocate it if needed
!
if (allocated(p%av)) then
if (size(p%av) < mld_bp_ilu_avsz_) then
do i=1, size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
endif
end if
endif
if (.not.allocated(p%av)) then
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_sp_get_nrows(a)
nztota = psb_sp_get_nnzeros(a)
if (present(blck)) then
nztota = nztota + psb_sp_get_nnzeros(blck)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ': out get_nnzeros',nztota,a%m,a%k,nrow_a
nrow_a = a%get_nrows()
nztota = a%get_nzeros()
if (present(blck)) then
nztota = nztota + blck%get_nzeros()
end if
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(mld_l_pr_)%m = n_row
p%av(mld_l_pr_)%k = n_row
p%av(mld_u_pr_)%m = n_row
p%av(mld_u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(mld_l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(mld_u_pr_),nztota,info)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call p%av(mld_l_pr_)%csall(n_row,n_row,info,nztota)
if (info == 0) call p%av(mld_u_pr_)%csall(n_row,n_row,info,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
endif
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
endif
select case(p%iprcparm(mld_sub_solve_))
case (mld_ilu_t_)
!
! ILU(k,t)
!
select case(p%iprcparm(mld_sub_solve_))
select case(p%iprcparm(mld_sub_fillin_))
case (mld_ilu_t_)
!
! ILU(k,t)
!
select case(p%iprcparm(mld_sub_fillin_))
case(:-1)
! Error: fill-in <= -1
call psb_errpush(30,name,i_err=(/3,p%iprcparm(mld_sub_fillin_),0,0,0/))
goto 9999
case(:-1)
! Error: fill-in <= -1
call psb_errpush(30,name,i_err=(/3,p%iprcparm(mld_sub_fillin_),0,0,0/))
goto 9999
case(0:)
! Fill-in >= 0
call mld_ilut_fact(p%iprcparm(mld_sub_fillin_),p%rprcparm(mld_sub_iluthrs_),&
& a, p%av(mld_l_pr_),p%av(mld_u_pr_),p%d,info,blck=blck)
end select
if(info/=0) then
info=4010
ch_err='mld_ilut_fact'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(0:)
! Fill-in >= 0
call mld_ilut_fact(p%iprcparm(mld_sub_fillin_),p%rprcparm(mld_sub_iluthrs_),&
& a, p%av(mld_l_pr_),p%av(mld_u_pr_),p%d,info,blck=blck)
end select
if(info/=0) then
info=4010
ch_err='mld_ilut_fact'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(mld_ilu_n_,mld_milu_n_)
!
! ILU(k) and MILU(k)
!
select case(p%iprcparm(mld_sub_fillin_))
case(:-1)
! Error: fill-in <= -1
call psb_errpush(30,name,i_err=(/3,p%iprcparm(mld_sub_fillin_),0,0,0/))
goto 9999
case(0)
! Fill-in 0
! Separate implementation of ILU(0) for better performance.
! There seems to be a problem with the separate implementation of MILU(0),
! contained into mld_ilu0_fact. This must be investigated. For the time being,
! resort to the implementation of MILU(k) with k=0.
if (p%iprcparm(mld_sub_solve_) == mld_ilu_n_) then
call mld_ilu0_fact(p%iprcparm(mld_sub_solve_),a,p%av(mld_l_pr_),p%av(mld_u_pr_),&
& p%d,info,blck=blck)
else
case(mld_ilu_n_,mld_milu_n_)
!
! ILU(k) and MILU(k)
!
select case(p%iprcparm(mld_sub_fillin_))
case(:-1)
! Error: fill-in <= -1
call psb_errpush(30,name,i_err=(/3,p%iprcparm(mld_sub_fillin_),0,0,0/))
goto 9999
case(0)
! Fill-in 0
! Separate implementation of ILU(0) for better performance.
! There seems to be a problem with the separate implementation of MILU(0),
! contained into mld_ilu0_fact. This must be investigated. For the time being,
! resort to the implementation of MILU(k) with k=0.
if (p%iprcparm(mld_sub_solve_) == mld_ilu_n_) then
call mld_ilu0_fact(p%iprcparm(mld_sub_solve_),a,p%av(mld_l_pr_),p%av(mld_u_pr_),&
& p%d,info,blck=blck,upd=upd)
else
call mld_iluk_fact(p%iprcparm(mld_sub_fillin_),p%iprcparm(mld_sub_solve_),&
& a,p%av(mld_l_pr_),p%av(mld_u_pr_),p%d,info,blck=blck)
endif
case(1:)
! Fill-in >= 1
! The same routine implements both ILU(k) and MILU(k)
call mld_iluk_fact(p%iprcparm(mld_sub_fillin_),p%iprcparm(mld_sub_solve_),&
& a,p%av(mld_l_pr_),p%av(mld_u_pr_),p%d,info,blck=blck)
endif
case(1:)
! Fill-in >= 1
! The same routine implements both ILU(k) and MILU(k)
call mld_iluk_fact(p%iprcparm(mld_sub_fillin_),p%iprcparm(mld_sub_solve_),&
& a,p%av(mld_l_pr_),p%av(mld_u_pr_),p%d,info,blck=blck)
end select
if (info/=0) then
info=4010
ch_err='mld_iluk_fact'
call psb_errpush(info,name,a_err=ch_err)
end select
if (info/=0) then
info=4010
ch_err='mld_iluk_fact'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case default
! If we end up here, something was wrong up in the call chain.
call psb_errpush(4000,name)
goto 9999
end if
case default
! If we end up here, something was wrong up in the call chain.
call psb_errpush(4000,name)
goto 9999
end select
else
! Here we should add checks for reuse of L and U.
! For the time being just throw an error.
info = 31
call psb_errpush(info, name, i_err=(/3,0,0,0,0/),a_err=upd)
goto 9999
end select
!
! What is an update of a factorization??
! A first attempt could be to reuse EXACTLY the existing indices
! as if it was an ILU(0) (since, effectively, the sparsity pattern
! should not grow beyond what is already there).
!
call mld_ilu0_fact(p%iprcparm(mld_sub_solve_),a,&
& p%av(mld_l_pr_),p%av(mld_u_pr_),&
& p%d,info,blck=blck,upd=upd)
if (psb_sp_getifld(psb_upd_,p%av(mld_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_u_pr_),info)
endif
end if
if (psb_sp_getifld(psb_upd_,p%av(mld_l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_l_pr_),info)
endif
call p%av(mld_l_pr_)%trim()
call p%av(mld_u_pr_)%trim()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

262
mlprec/mld_diluk_fact.f90
Unescape Escape View File

@ -106,14 +106,15 @@ subroutine mld_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
! Arguments
integer, intent(in) :: fill_in, ialg
integer, intent(out) :: info
type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u
type(psb_dspmat_type),intent(in), optional, target :: blck
type(psb_d_sparse_mat),intent(in) :: a
type(psb_d_sparse_mat),intent(inout) :: l,u
type(psb_d_sparse_mat),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:)
! Local Variables
integer :: l1, l2, m, err_act
type(psb_dspmat_type), pointer :: blck_
type(psb_d_sparse_mat), pointer :: blck_
type(psb_d_csr_sparse_mat) :: ll, uu
character(len=20) :: name, ch_err
name='mld_diluk_fact'
@ -127,26 +128,32 @@ subroutine mld_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
blck_ => blck
else
allocate(blck_,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
call psb_sp_all(0,0,blck_,1,info)
if (info == 0) call blck_%csall(0,0,info,1)
if (info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='csall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
endif
m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
& (m > size(d)) ) then
write(0,*) 'Wrong allocation status for L,D,U? ',&
& l%get_nrows(),size(d),u%get_nrows()
info = -1
return
end if
call l%mv_to(ll)
call u%mv_to(uu)
!
! Compute the ILU(k) or the MILU(k) factorization, depending on ialg
!
call mld_diluk_factint(fill_in,ialg,m,a,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
call mld_diluk_factint(fill_in,ialg,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info)
if (info /= 0) then
info=4010
ch_err='mld_diluk_factint'
@ -157,33 +164,32 @@ subroutine mld_diluk_fact(fill_in,ialg,a,l,u,d,info,blck)
!
! Store information on the L and U sparse matrices
!
l%infoa(1) = l1
l%fida = 'CSR'
l%descra = 'TLU'
u%infoa(1) = l2
u%fida = 'CSR'
u%descra = 'TUU'
l%m = m
l%k = m
u%m = m
u%k = m
call l%mv_from(ll)
call l%set_triangle()
call l%set_unit()
call l%set_lower()
call u%mv_from(uu)
call u%set_triangle()
call u%set_unit()
call u%set_upper()
!
! Nullify the pointer / deallocate the memory
! Nullify pointer / deallocate memory
!
if (present(blck)) then
blck_ => null()
else
call psb_sp_free(blck_,info)
if (info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
call blck_%free()
deallocate(blck_,stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
deallocate(blck_)
endif
call psb_erractionrestore(err_act)
return
@ -248,43 +254,43 @@ contains
! lia2 - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the L factor in laspk, according to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! uval - real(psb_dpk_), dimension(:), input/output.
! The U factor in the incomplete factorization.
! The entries of U are stored according to the CSR format.
! uia1 - integer, dimension(:), input/output.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the U factor,
! according to the CSR storage format.
! uia2 - integer, dimension(:), input/output.
! uirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the U factor in uaspk, according to the CSR storage format.
! of the U factor in uval, according to the CSR storage format.
! l1 - integer, output
! The number of nonzero entries in laspk.
! l2 - integer, output
! The number of nonzero entries in uaspk.
! The number of nonzero entries in uval.
! info - integer, output.
! Error code.
!
subroutine mld_diluk_factint(fill_in,ialg,m,a,b,&
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
subroutine mld_diluk_factint(fill_in,ialg,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info)
use psb_base_mod
implicit none
! Arguments
integer, intent(in) :: fill_in, ialg
type(psb_dspmat_type), intent(in) :: a,b
integer, intent(inout) :: m,l1,l2,info
integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:)
real(psb_dpk_), allocatable, intent(inout) :: laspk(:),uaspk(:)
integer, intent(in) :: fill_in, ialg
type(psb_d_sparse_mat),intent(in) :: a,b
integer,intent(inout) :: l1,l2,info
integer, allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_dpk_), intent(inout) :: d(:)
! Local variables
integer :: ma,mb,i, ktrw,err_act,nidx
integer :: ma,mb,i, ktrw,err_act,nidx, m
integer, allocatable :: uplevs(:), rowlevs(:),idxs(:)
real(psb_dpk_), allocatable :: row(:)
real(psb_dpk_), allocatable :: row(:)
type(psb_int_heap) :: heap
type(psb_dspmat_type) :: trw
type(psb_d_coo_sparse_mat) :: trw
character(len=20), parameter :: name='mld_diluk_factint'
character(len=20) :: ch_err
@ -292,6 +298,7 @@ contains
info=0
call psb_erractionsave(err_act)
select case(ialg)
case(mld_ilu_n_,mld_milu_n_)
! Ok
@ -306,16 +313,17 @@ contains
goto 9999
end if
ma = a%m
mb = b%m
ma = a%get_nrows()
mb = b%get_nrows()
m = ma+mb
!
! Allocate a temporary buffer for the iluk_copyin function
!
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
call trw%allocate(0,0,1)
if (info==0) call psb_ensure_size(m+1,lirp,info)
if (info==0) call psb_ensure_size(m+1,uirp,info)
if (info /= 0) then
info=4010
@ -325,14 +333,14 @@ contains
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
lirp(1) = 1
uirp(1) = 1
!
! Allocate memory to hold the entries of a row and the corresponding
! fill levels
!
allocate(uplevs(size(uaspk)),rowlevs(m),row(m),stat=info)
allocate(uplevs(size(uval)),rowlevs(m),row(m),stat=info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='Allocate')
@ -375,12 +383,12 @@ contains
! do not have a lowlevs variable.
!
if (info == 0) call iluk_fact(fill_in,i,row,rowlevs,heap,&
& d,uia1,uia2,uaspk,uplevs,nidx,idxs,info)
& d,uja,uirp,uval,uplevs,nidx,idxs,info)
!
! Copy the row into laspk/d(i)/uaspk
! Copy the row into lval/d(i)/uval
!
if (info == 0) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,&
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info)
& l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info)
if (info /= 0) then
info=4001
call psb_errpush(info,name,a_err='Copy/factor loop')
@ -397,7 +405,7 @@ contains
call psb_errpush(info,name,a_err='Deallocate')
goto 9999
end if
if (info == 0) call psb_sp_free(trw,info)
if (info == 0) call trw%free()
if (info /= 0) then
info=4010
ch_err='psb_sp_free'
@ -488,8 +496,8 @@ contains
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(inout) :: trw
type(psb_d_sparse_mat), intent(in) :: a
type(psb_d_coo_sparse_mat), intent(inout) :: trw
integer, intent(in) :: i,m,jmin,jmax
integer, intent(inout) :: ktrw,info
integer, intent(inout) :: rowlevs(:)
@ -497,8 +505,8 @@ contains
type(psb_int_heap), intent(inout) :: heap
! Local variables
integer :: k,j,irb,err_act
integer, parameter :: nrb=16
integer :: k,j,irb,err_act,nz
integer, parameter :: nrb=40
character(len=20), parameter :: name='iluk_copyin'
character(len=20) :: ch_err
@ -507,22 +515,22 @@ contains
call psb_erractionsave(err_act)
call psb_init_heap(heap,info)
if (psb_toupper(a%fida)=='CSR') then
select type (aa=> a%a)
type is (psb_d_csr_sparse_mat)
!
! Take a fast shortcut if the matrix is stored in CSR format
!
do j = a%ia2(i), a%ia2(i+1) - 1
k = a%ia1(j)
do j = aa%irp(i), aa%irp(i+1) - 1
k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then
row(k) = a%aspk(j)
row(k) = aa%val(j)
rowlevs(k) = 0
call psb_insert_heap(k,heap,info)
end if
end do
else
class default
!
! Otherwise use psb_sp_getblk, slower but able (in principle) of
@ -534,7 +542,7 @@ contains
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
call aa%csget(i,i+irb-1,trw,info)
if (info /= 0) then
info=4010
ch_err='psb_sp_getblk'
@ -543,19 +551,19 @@ contains
end if
ktrw=1
end if
nz = trw%get_nzeros()
do
if (ktrw > trw%infoa(psb_nnz_)) exit
if (trw%ia1(ktrw) > i) exit
k = trw%ia2(ktrw)
if (ktrw > nz) exit
if (trw%ia(ktrw) > i) exit
k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%aspk(ktrw)
row(k) = trw%val(ktrw)
rowlevs(k) = 0
call psb_insert_heap(k,heap,info)
end if
ktrw = ktrw + 1
enddo
end if
end select
call psb_erractionrestore(err_act)
return
@ -611,17 +619,17 @@ contains
! d - real(psb_dpk_), input.
! The inverse of the diagonal entries of the part of the U factor
! above the current row (see iluk_copyout).
! uia1 - integer, dimension(:), input.
! uja - integer, dimension(:), input.
! The column indices of the nonzero entries of the part of the U
! factor above the current row, stored in uaspk row by row (see
! factor above the current row, stored in uval row by row (see
! iluk_copyout, called by mld_diluk_factint), according to the CSR
! storage format.
! uia2 - integer, dimension(:), input.
! uirp - integer, dimension(:), input.
! The indices identifying the first nonzero entry of each row of
! the U factor above the current row, stored in uaspk row by row
! the U factor above the current row, stored in uval row by row
! (see iluk_copyout, called by mld_diluk_factint), according to
! the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input.
! uval - real(psb_dpk_), dimension(:), input.
! The entries of the U factor above the current row (except the
! diagonal ones), stored according to the CSR format.
! uplevs - integer, dimension(:), input.
@ -638,7 +646,7 @@ contains
! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine.
!
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uia1,uia2,uaspk,uplevs,nidx,idxs,info)
subroutine iluk_fact(fill_in,i,row,rowlevs,heap,d,uja,uirp,uval,uplevs,nidx,idxs,info)
use psb_base_mod
@ -650,8 +658,8 @@ contains
integer, intent(inout) :: nidx,info
integer, intent(inout) :: rowlevs(:)
integer, allocatable, intent(inout) :: idxs(:)
integer, intent(inout) :: uia1(:),uia2(:),uplevs(:)
real(psb_dpk_), intent(inout) :: row(:), uaspk(:),d(:)
integer, intent(inout) :: uja(:),uirp(:),uplevs(:)
real(psb_dpk_), intent(inout) :: row(:), uval(:),d(:)
! Local variables
integer :: k,j,lrwk,jj,lastk, iret
@ -695,8 +703,8 @@ contains
row(k) = row(k) * d(k) ! d(k) == 1/a(k,k)
lrwk = rowlevs(k)
do jj=uia2(k),uia2(k+1)-1
j = uia1(jj)
do jj=uirp(k),uirp(k+1)-1
j = uja(jj)
if (j<=k) then
info = -i
return
@ -716,7 +724,7 @@ contains
!
! Update row(j) and the corresponding fill level
!
row(j) = row(j) - rwk * uaspk(jj)
row(j) = row(j) - rwk * uval(jj)
rowlevs(j) = min(rowlevs(j),lrwk+uplevs(jj)+1)
end do
@ -731,19 +739,19 @@ contains
! Note: internal subroutine of mld_diluk_fact
!
! This routine copies a matrix row, computed by iluk_fact by applying an
! elimination step of the ILU(k) factorization, into the arrays laspk, uaspk,
! elimination step of the ILU(k) factorization, into the arrays lval, uval,
! d, corresponding to the L factor, the U factor and the diagonal of U,
! respectively.
!
! Note that
! - the part of the row stored into uaspk is scaled by the corresponding diagonal
! - the part of the row stored into uval is scaled by the corresponding diagonal
! entry, according to the LDU form of the incomplete factorization;
! - the inverse of the diagonal entries of U is actually stored into d; this is
! then managed in the solve stage associated to the ILU(k)/MILU(k) factorization;
! - if the MILU(k) factorization has been required (ialg == mld_milu_n_), the
! row entries discarded because their fill levels are too high are added to
! the diagonal entry of the row;
! - the row entries are stored in laspk and uaspk according to the CSR format;
! - the row entries are stored in lval and uval according to the CSR format;
! - the arrays row and rowlevs are re-initialized for future use in mld_iluk_fact
! (see also iluk_copyin and iluk_fact).
!
@ -781,32 +789,32 @@ contains
! examined during the elimination step carried out by the routine
! iluk_fact.
! l1 - integer, input/output.
! Pointer to the last occupied entry of laspk.
! Pointer to the last occupied entry of lval.
! l2 - integer, input/output.
! Pointer to the last occupied entry of uaspk.
! lia1 - integer, dimension(:), input/output.
! Pointer to the last occupied entry of uval.
! lja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the L factor,
! copied in laspk row by row (see mld_diluk_factint), according
! copied in lval row by row (see mld_diluk_factint), according
! to the CSR storage format.
! lia2 - integer, dimension(:), input/output.
! lirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the L factor, copied in laspk row by row (see
! of the L factor, copied in lval row by row (see
! mld_diluk_factint), according to the CSR storage format.
! laspk - real(psb_dpk_), dimension(:), input/output.
! lval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! L factor are copied.
! d - real(psb_dpk_), dimension(:), input/output.
! The array where the inverse of the diagonal entry of the
! row is copied (only d(i) is used by the routine).
! uia1 - integer, dimension(:), input/output.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the U factor
! copied in uaspk row by row (see mld_diluk_factint), according
! copied in uval row by row (see mld_diluk_factint), according
! to the CSR storage format.
! uia2 - integer, dimension(:), input/output.
! uirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the U factor copied in uaspk row by row (see
! of the U factor copied in uval row by row (see
! mld_dilu_fctint), according to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! uval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! U factor are copied.
! uplevs - integer, dimension(:), input.
@ -814,18 +822,18 @@ contains
! U factor above the current row.
!
subroutine iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,&
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,uplevs,info)
& l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info)
use psb_base_mod
implicit none
! Arguments
integer, intent(in) :: fill_in, ialg, i, m, nidx
integer, intent(inout) :: l1, l2, info
integer, intent(inout) :: rowlevs(:), idxs(:)
integer, allocatable, intent(inout) :: uia1(:), uia2(:), lia1(:), lia2(:),uplevs(:)
real(psb_dpk_), allocatable, intent(inout) :: uaspk(:), laspk(:)
integer, intent(in) :: fill_in, ialg, i, m, nidx
integer, intent(inout) :: l1, l2, info
integer, intent(inout) :: rowlevs(:), idxs(:)
integer, allocatable, intent(inout) :: uja(:), uirp(:), lja(:), lirp(:),uplevs(:)
real(psb_dpk_), allocatable, intent(inout) :: uval(:), lval(:)
real(psb_dpk_), intent(inout) :: row(:), d(:)
! Local variables
@ -849,21 +857,21 @@ contains
!
if (rowlevs(j) <= fill_in) then
l1 = l1 + 1
if (size(laspk) < l1) then
if (size(lval) < l1) then
!
! Figure out a good reallocation size!
!
isz = (max((l1/i)*m,int(1.2*l1),l1+100))
call psb_realloc(isz,laspk,info)
if (info == 0) call psb_realloc(isz,lia1,info)
call psb_realloc(isz,lval,info)
if (info == 0) call psb_realloc(isz,lja,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
end if
lia1(l1) = j
laspk(l1) = row(j)
lja(l1) = j
lval(l1) = row(j)
else if (ialg == mld_milu_n_) then
!
! MILU(k): add discarded entries to the diagonal one
@ -891,13 +899,13 @@ contains
!
if (rowlevs(j) <= fill_in) then
l2 = l2 + 1
if (size(uaspk) < l2) then
if (size(uval) < l2) then
!
! Figure out a good reallocation size!
!
isz = max((l2/i)*m,int(1.2*l2),l2+100)
call psb_realloc(isz,uaspk,info)
if (info == 0) call psb_realloc(isz,uia1,info)
call psb_realloc(isz,uval,info)
if (info == 0) call psb_realloc(isz,uja,info)
if (info == 0) call psb_realloc(isz,uplevs,info,pad=(m+1))
if (info /= 0) then
info=4010
@ -905,8 +913,8 @@ contains
goto 9999
end if
end if
uia1(l2) = j
uaspk(l2) = row(j)
uja(l2) = j
uval(l2) = row(j)
uplevs(l2) = rowlevs(j)
else if (ialg == mld_milu_n_) then
!
@ -924,10 +932,10 @@ contains
!
! Store the pointers to the first non occupied entry of in
! laspk and uaspk
! lval and uval
!
lia2(i+1) = l1 + 1
uia2(i+1) = l2 + 1
lirp(i+1) = l1 + 1
uirp(i+1) = l2 + 1
!
! Check the pivot size
@ -951,8 +959,8 @@ contains
!
! Scale the upper part
!
do j=uia2(i), uia2(i+1)-1
uaspk(j) = d(i)*uaspk(j)
do j=uirp(i), uirp(i+1)-1
uval(j) = d(i)*uval(j)
end do
call psb_erractionrestore(err_act)

297
mlprec/mld_dilut_fact.f90
Unescape Escape View File

@ -103,15 +103,15 @@ subroutine mld_dilut_fact(fill_in,thres,a,l,u,d,info,blck)
integer, intent(in) :: fill_in
real(psb_dpk_), intent(in) :: thres
integer, intent(out) :: info
type(psb_dspmat_type),intent(in) :: a
type(psb_dspmat_type),intent(inout) :: l,u
real(psb_dpk_), intent(inout) :: d(:)
type(psb_dspmat_type),intent(in), optional, target :: blck
type(psb_d_sparse_mat),intent(in) :: a
type(psb_d_sparse_mat),intent(inout) :: l,u
type(psb_d_sparse_mat),intent(in), optional, target :: blck
real(psb_dpk_), intent(inout) :: d(:)
! Local Variables
integer :: l1, l2, m, err_act
type(psb_dspmat_type), pointer :: blck_
type(psb_d_sparse_mat), pointer :: blck_
type(psb_d_csr_sparse_mat) :: ll, uu
character(len=20) :: name, ch_err
name='mld_dilut_fact'
@ -130,26 +130,32 @@ subroutine mld_dilut_fact(fill_in,thres,a,l,u,d,info,blck)
blck_ => blck
else
allocate(blck_,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
call psb_sp_all(0,0,blck_,1,info)
if (info == 0) call blck_%csall(0,0,info,1)
if (info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
info=4010
ch_err='csall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
endif
m = a%get_nrows() + blck_%get_nrows()
if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.&
& (m > size(d)) ) then
write(0,*) 'Wrong allocation status for L,D,U? ',&
& l%get_nrows(),size(d),u%get_nrows()
info = -1
return
end if
call l%mv_to(ll)
call u%mv_to(uu)
!
! Compute the ILU(k,t) factorization
!
call mld_dilut_factint(fill_in,thres,m,a,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
call mld_dilut_factint(fill_in,thres,a,blck_,&
& d,ll%val,ll%ja,ll%irp,uu%val,uu%ja,uu%irp,l1,l2,info)
if (info /= 0) then
info=4010
ch_err='mld_dilut_factint'
@ -160,31 +166,29 @@ subroutine mld_dilut_fact(fill_in,thres,a,l,u,d,info,blck)
!
! Store information on the L and U sparse matrices
!
l%infoa(1) = l1
l%fida = 'CSR'
l%descra = 'TLU'
u%infoa(1) = l2
u%fida = 'CSR'
u%descra = 'TUU'
l%m = m
l%k = m
u%m = m
u%k = m
call l%mv_from(ll)
call l%set_triangle()
call l%set_unit()
call l%set_lower()
call u%mv_from(uu)
call u%set_triangle()
call u%set_unit()
call u%set_upper()
!
! Nullify the pointer / deallocate the memory
! Nullify pointer / deallocate memory
!
if (present(blck)) then
blck_ => null()
else
call psb_sp_free(blck_,info)
if (info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
call blck_%free()
deallocate(blck_,stat=info)
if(info.ne.0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
deallocate(blck_)
endif
call psb_erractionrestore(err_act)
@ -241,32 +245,32 @@ contains
! d - real(psb_dpk_), dimension(:), output.
! The inverse of the diagonal entries of the U factor in the incomplete
! factorization.
! laspk - real(psb_dpk_), dimension(:), input/output.
! lval - real(psb_dpk_), dimension(:), input/output.
! The L factor in the incomplete factorization.
! lia1 - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the L factor,
! according to the CSR storage format.
! lia2 - integer, dimension(:), input/output.
! lirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the L factor in laspk, according to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! of the L factor in lval, according to the CSR storage format.
! uval - real(psb_dpk_), dimension(:), input/output.
! The U factor in the incomplete factorization.
! The entries of U are stored according to the CSR format.
! uia1 - integer, dimension(:), input/output.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the U factor,
! according to the CSR storage format.
! uia2 - integer, dimension(:), input/output.
! uirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the U factor in uaspk, according to the CSR storage format.
! of the U factor in uval, according to the CSR storage format.
! l1 - integer, output
! The number of nonzero entries in laspk.
! The number of nonzero entries in lval.
! l2 - integer, output
! The number of nonzero entries in uaspk.
! The number of nonzero entries in uval.
! info - integer, output.
! Error code.
!
subroutine mld_dilut_factint(fill_in,thres,m,a,b,&
& d,laspk,lia1,lia2,uaspk,uia1,uia2,l1,l2,info)
subroutine mld_dilut_factint(fill_in,thres,a,b,&
& d,lval,lja,lirp,uval,uja,uirp,l1,l2,info)
use psb_base_mod
@ -275,19 +279,19 @@ contains
! Arguments
integer, intent(in) :: fill_in
real(psb_dpk_), intent(in) :: thres
type(psb_dspmat_type), intent(in) :: a,b
integer, intent(inout) :: m,l1,l2,info
integer, allocatable, intent(inout) :: lia1(:),lia2(:),uia1(:),uia2(:)
real(psb_dpk_), allocatable, intent(inout) :: laspk(:),uaspk(:)
type(psb_d_sparse_mat),intent(in) :: a,b
integer,intent(inout) :: l1,l2,info
integer, allocatable, intent(inout) :: lja(:),lirp(:),uja(:),uirp(:)
real(psb_dpk_), allocatable, intent(inout) :: lval(:),uval(:)
real(psb_dpk_), intent(inout) :: d(:)
! Local Variables
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb
integer :: i, ktrw,err_act,nidx,nlw,nup,jmaxup, ma, mb, m
real(psb_dpk_) :: nrmi
integer, allocatable :: idxs(:)
real(psb_dpk_), allocatable :: row(:)
type(psb_int_heap) :: heap
type(psb_dspmat_type) :: trw
type(psb_d_coo_sparse_mat) :: trw
character(len=20), parameter :: name='mld_dilut_factint'
character(len=20) :: ch_err
@ -296,16 +300,16 @@ contains
call psb_erractionsave(err_act)
ma = a%m
mb = b%m
ma = a%get_nrows()
mb = b%get_nrows()
m = ma+mb
!
! Allocate a temporary buffer for the ilut_copyin function
!
call psb_sp_all(0,0,trw,1,info)
if (info==0) call psb_ensure_size(m+1,lia2,info)
if (info==0) call psb_ensure_size(m+1,uia2,info)
call trw%allocate(0,0,1)
if (info==0) call psb_ensure_size(m+1,lirp,info)
if (info==0) call psb_ensure_size(m+1,uirp,info)
if (info /= 0) then
info=4010
@ -315,8 +319,8 @@ contains
l1=0
l2=0
lia2(1) = 1
uia2(1) = 1
lirp(1) = 1
uirp(1) = 1
!
! Allocate memory to hold the entries of a row
@ -354,12 +358,12 @@ contains
! Do an elimination step on current row
!
if (info == 0) call ilut_fact(thres,i,nrmi,row,heap,&
& d,uia1,uia2,uaspk,nidx,idxs,info)
& d,uja,uirp,uval,nidx,idxs,info)
!
! Copy the row into laspk/d(i)/uaspk
! Copy the row into lval/d(i)/uval
!
if (info == 0) call ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row,nidx,idxs,&
& l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
& l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
if (info /= 0) then
info=4001
@ -378,7 +382,7 @@ contains
call psb_errpush(info,name,a_err='Deallocate')
goto 9999
end if
if (info == 0) call psb_sp_free(trw,info)
if (info == 0) call trw%free()
if (info /= 0) then
info=4010
ch_err='psb_sp_free'
@ -482,16 +486,16 @@ contains
subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,nrmi,row,heap,ktrw,trw,info)
use psb_base_mod
implicit none
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(inout) :: trw
integer, intent(in) :: i, m,jmin,jmax,jd
integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
real(psb_dpk_), intent(inout) :: nrmi,row(:)
type(psb_int_heap), intent(inout) :: heap
type(psb_d_sparse_mat), intent(in) :: a
type(psb_d_coo_sparse_mat), intent(inout) :: trw
integer, intent(in) :: i, m,jmin,jmax,jd
integer, intent(inout) :: ktrw,nlw,nup,jmaxup,info
real(psb_dpk_), intent(inout) :: nrmi,row(:)
type(psb_int_heap), intent(inout) :: heap
integer :: k,j,irb,kin,nz
integer, parameter :: nrb=16
real(psb_dpk_) :: dmaxup
integer, parameter :: nrb=40
real(psb_dpk_) :: dmaxup
real(psb_dpk_), external :: dnrm2
character(len=20), parameter :: name='mld_dilut_factint'
@ -518,23 +522,19 @@ contains
jmaxup = 0
dmaxup = dzero
nrmi = dzero
if (psb_toupper(a%fida)=='CSR') then
select type (aa=> a%a)
type is (psb_d_csr_sparse_mat)
!
! Take a fast shortcut if the matrix is stored in CSR format
!
do j = a%ia2(i), a%ia2(i+1) - 1
k = a%ia1(j)
!
do j = aa%irp(i), aa%irp(i+1) - 1
k = aa%ja(j)
if ((jmin<=k).and.(k<=jmax)) then
row(k) = a%aspk(j)
row(k) = aa%val(j)
call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
if (info /= 0) exit
end if
if (k<jd) nlw = nlw + 1
if (k>jd) then
@ -545,9 +545,17 @@ contains
end if
end if
end do
nz = a%ia2(i+1) - a%ia2(i)
nrmi = dnrm2(nz,a%aspk(a%ia2(i)),ione)
else
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
nz = aa%irp(i+1) - aa%irp(i)
nrmi = dnrm2(nz,aa%val(aa%irp(i)),ione)
class default
!
! Otherwise use psb_sp_getblk, slower but able (in principle) of
@ -559,7 +567,7 @@ contains
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
call aa%csget(i,i+irb-1,trw,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_sp_getblk')
@ -569,18 +577,16 @@ contains
end if
kin = ktrw
nz = trw%get_nzeros()
do
if (ktrw > trw%infoa(psb_nnz_)) exit
if (trw%ia1(ktrw) > i) exit
k = trw%ia2(ktrw)
if (ktrw > nz) exit
if (trw%ia(ktrw) > i) exit
k = trw%ja(ktrw)
if ((jmin<=k).and.(k<=jmax)) then
row(k) = trw%aspk(ktrw)
row(k) = trw%val(ktrw)
call psb_insert_heap(k,heap,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_insert_heap')
goto 9999
end if
if (info /= 0) exit
end if
if (k<jd) nlw = nlw + 1
if (k>jd) then
@ -593,8 +599,9 @@ contains
ktrw = ktrw + 1
enddo
nz = ktrw - kin
nrmi = dnrm2(nz,trw%aspk(kin),ione)
end if
nrmi = dnrm2(nz,trw%val(kin),ione)
end select
call psb_erractionrestore(err_act)
return
@ -644,17 +651,17 @@ contains
! d - real(psb_dpk_), input.
! The inverse of the diagonal entries of the part of the U factor
! above the current row (see ilut_copyout).
! uia1 - integer, dimension(:), input.
! uja - integer, dimension(:), input.
! The column indices of the nonzero entries of the part of the U
! factor above the current row, stored in uaspk row by row (see
! factor above the current row, stored in uval row by row (see
! ilut_copyout, called by mld_dilut_factint), according to the CSR
! storage format.
! uia2 - integer, dimension(:), input.
! uirp - integer, dimension(:), input.
! The indices identifying the first nonzero entry of each row of
! the U factor above the current row, stored in uaspk row by row
! the U factor above the current row, stored in uval row by row
! (see ilut_copyout, called by mld_dilut_factint), according to
! the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input.
! uval - real(psb_dpk_), dimension(:), input.
! The entries of the U factor above the current row (except the
! diagonal ones), stored according to the CSR format.
! nidx - integer, output.
@ -668,7 +675,7 @@ contains
! Note: this argument is intent(inout) and not only intent(out)
! to retain its allocation, done by this routine.
!
subroutine ilut_fact(thres,i,nrmi,row,heap,d,uia1,uia2,uaspk,nidx,idxs,info)
subroutine ilut_fact(thres,i,nrmi,row,heap,d,uja,uirp,uval,nidx,idxs,info)
use psb_base_mod
@ -680,8 +687,8 @@ contains
integer, intent(inout) :: nidx,info
real(psb_dpk_), intent(in) :: thres,nrmi
integer, allocatable, intent(inout) :: idxs(:)
integer, intent(inout) :: uia1(:),uia2(:)
real(psb_dpk_), intent(inout) :: row(:), uaspk(:),d(:)
integer, intent(inout) :: uja(:),uirp(:)
real(psb_dpk_), intent(inout) :: row(:), uval(:),d(:)
! Local Variables
integer :: k,j,jj,lastk,iret
@ -725,8 +732,8 @@ contains
! Note: since U is scaled while copying it out (see ilut_copyout),
! we can use rwk in the update below.
!
do jj=uia2(k),uia2(k+1)-1
j = uia1(jj)
do jj=uirp(k),uirp(k+1)-1
j = uja(jj)
if (j<=k) then
info = -i
return
@ -735,7 +742,7 @@ contains
! Update row(j) and, if it is not to be discarded, insert
! its index into the heap for further processing.
!
row(j) = row(j) - rwk * uaspk(jj)
row(j) = row(j) - rwk * uval(jj)
if (abs(row(j)) < thres*nrmi) then
!
! Drop the entry.
@ -770,8 +777,8 @@ contains
! Note: internal subroutine of mld_dilut_fact
!
! This routine copies a matrix row, computed by ilut_fact by applying an
! elimination step of the ILU(k,t) factorization, into the arrays laspk,
! uaspk, d, corresponding to the L factor, the U factor and the diagonal
! elimination step of the ILU(k,t) factorization, into the arrays lval,
! uval, d, corresponding to the L factor, the U factor and the diagonal
! of U, respectively.
!
! Note that
@ -780,11 +787,11 @@ contains
! the 'lower part' of the row, and the nup+k ones in the 'upper part';
! - the entry in the upper part of the row which has maximum absolute value
! in the original matrix is included in the above nup+k entries anyway;
! - the part of the row stored into uaspk is scaled by the corresponding
! - the part of the row stored into uval is scaled by the corresponding
! diagonal entry, according to the LDU form of the incomplete factorization;
! - the inverse of the diagonal entries of U is actually stored into d; this
! is then managed in the solve stage associated to the ILU(k,t) factorization;
! - the row entries are stored in laspk and uaspk according to the CSR format;
! - the row entries are stored in lval and uval according to the CSR format;
! - the array row is re-initialized for future use in mld_ilut_fact(see also
! ilut_copyin and ilut_fact).
!
@ -824,37 +831,37 @@ contains
! examined during the elimination step carried out by the routine
! ilut_fact.
! l1 - integer, input/output.
! Pointer to the last occupied entry of laspk.
! Pointer to the last occupied entry of lval.
! l2 - integer, input/output.
! Pointer to the last occupied entry of uaspk.
! lia1 - integer, dimension(:), input/output.
! Pointer to the last occupied entry of uval.
! lja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the L factor,
! copied in laspk row by row (see mld_dilut_factint), according
! copied in lval row by row (see mld_dilut_factint), according
! to the CSR storage format.
! lia2 - integer, dimension(:), input/output.
! lirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the L factor, copied in laspk row by row (see
! of the L factor, copied in lval row by row (see
! mld_dilut_factint), according to the CSR storage format.
! laspk - real(psb_dpk_), dimension(:), input/output.
! lval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! L factor are copied.
! d - real(psb_dpk_), dimension(:), input/output.
! The array where the inverse of the diagonal entry of the
! row is copied (only d(i) is used by the routine).
! uia1 - integer, dimension(:), input/output.
! uja - integer, dimension(:), input/output.
! The column indices of the nonzero entries of the U factor
! copied in uaspk row by row (see mld_dilut_factint), according
! copied in uval row by row (see mld_dilut_factint), according
! to the CSR storage format.
! uia2 - integer, dimension(:), input/output.
! uirp - integer, dimension(:), input/output.
! The indices identifying the first nonzero entry of each row
! of the U factor copied in uaspk row by row (see
! of the U factor copied in uval row by row (see
! mld_dilu_fctint), according to the CSR storage format.
! uaspk - real(psb_dpk_), dimension(:), input/output.
! uval - real(psb_dpk_), dimension(:), input/output.
! The array where the entries of the row corresponding to the
! U factor are copied.
!
subroutine ilut_copyout(fill_in,thres,i,m,nlw,nup,jmaxup,nrmi,row, &
& nidx,idxs,l1,l2,lia1,lia2,laspk,d,uia1,uia2,uaspk,info)
& nidx,idxs,l1,l2,lja,lirp,lval,d,uja,uirp,uval,info)
use psb_base_mod
@ -864,9 +871,9 @@ contains
integer, intent(in) :: fill_in,i,m,nidx,nlw,nup,jmaxup
integer, intent(in) :: idxs(:)
integer, intent(inout) :: l1,l2, info
integer, allocatable, intent(inout) :: uia1(:),uia2(:), lia1(:),lia2(:)
integer, allocatable, intent(inout) :: uja(:),uirp(:), lja(:),lirp(:)
real(psb_dpk_), intent(in) :: thres,nrmi
real(psb_dpk_),allocatable, intent(inout) :: uaspk(:), laspk(:)
real(psb_dpk_),allocatable, intent(inout) :: uval(:), lval(:)
real(psb_dpk_), intent(inout) :: row(:), d(:)
! Local variables
@ -965,21 +972,21 @@ contains
!
do k=1,nz
l1 = l1 + 1
if (size(laspk) < l1) then
if (size(lval) < l1) then
!
! Figure out a good reallocation size!
!
isz = (max((l1/i)*m,int(1.2*l1),l1+100))
call psb_realloc(isz,laspk,info)
if (info == 0) call psb_realloc(isz,lia1,info)
call psb_realloc(isz,lval,info)
if (info == 0) call psb_realloc(isz,lja,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
end if
lia1(l1) = xwid(k)
laspk(l1) = xw(indx(k))
lja(l1) = xwid(k)
lval(l1) = xw(indx(k))
end do
!
@ -1111,21 +1118,21 @@ contains
!
do k=1,nz
l2 = l2 + 1
if (size(uaspk) < l2) then
if (size(uval) < l2) then
!
! Figure out a good reallocation size!
!
isz = max((l2/i)*m,int(1.2*l2),l2+100)
call psb_realloc(isz,uaspk,info)
if (info == 0) call psb_realloc(isz,uia1,info)
call psb_realloc(isz,uval,info)
if (info == 0) call psb_realloc(isz,uja,info)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
end if
uia1(l2) = xwid(k)
uaspk(l2) = d(i)*xw(indx(k))
uja(l2) = xwid(k)
uval(l2) = d(i)*xw(indx(k))
end do
!
@ -1137,10 +1144,10 @@ contains
!
! Store the pointers to the first non occupied entry of in
! laspk and uaspk
! lval and uval
!
lia2(i+1) = l1 + 1
uia2(i+1) = l2 + 1
lirp(i+1) = l1 + 1
uirp(i+1) = l2 + 1
call psb_erractionrestore(err_act)
return

14
mlprec/mld_dmlprec_bld.f90
Unescape Escape View File

@ -73,10 +73,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
Implicit None
! Arguments
type(psb_dspmat_type),intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type),intent(inout),target :: p
integer, intent(out) :: info
type(psb_d_sparse_mat),intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type),intent(inout),target :: p
integer, intent(out) :: info
!!$ character, intent(in), optional :: upd
! Local Variables
@ -346,9 +346,9 @@ contains
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) return
end if
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do
!!$ do k=1,size(p%av)
!!$ call psb_nullify_sp(p%av(k))
!!$ end do
end subroutine init_baseprec_av

4
mlprec/mld_dprecaply.f90
Unescape Escape View File

@ -74,7 +74,7 @@
subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
use mld_inner_mod
use mld_inner_mod, mld_protect_name => mld_dprecaply
use mld_prec_mod, mld_protect_name => mld_dprecaply
implicit none
@ -206,7 +206,7 @@ end subroutine mld_dprecaply
subroutine mld_dprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod
use mld_inner_mod
use mld_inner_mod, mld_protect_name => mld_dprecaply1
use mld_prec_mod, mld_protect_name => mld_dprecaply1
implicit none

8
mlprec/mld_dprecbld.f90
Unescape Escape View File

@ -67,7 +67,7 @@ subroutine mld_dprecbld(a,desc_a,p,info)
Implicit None
! Arguments
type(psb_dspmat_type),intent(in), target :: a
type(psb_d_sparse_mat),intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type),intent(inout),target :: p
integer, intent(out) :: info
@ -231,9 +231,9 @@ contains
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) return
end if
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do
!!$ do k=1,size(p%av)
!!$ call psb_nullify_sp(p%av(k))
!!$ end do
end subroutine init_baseprec_av

39
mlprec/mld_dslu_bld.f90
Unescape Escape View File

@ -77,14 +77,14 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_d_sparse_mat), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(mld_dbaseprec_type), intent(inout) :: p
integer, intent(out) :: info
integer, intent(out) :: info
! Local variables
integer :: nzt,ictxt,me,np,err_act
character(len=20) :: name, ch_err
integer :: ictxt,me,np,err_act
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
info=0
@ -94,25 +94,28 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
if (psb_toupper(a%fida) /= 'CSR') then
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
!
! Compute the LU factorization
!
call mld_dslu_fact(a%m,nzt,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slu_ptr_),info)
select type(aa=>a%a)
type is (psb_d_csr_sparse_mat)
call mld_dslu_fact(aa%get_nrows(),aa%get_nzeros(),&
& aa%val,aa%ja,aa%irp,p%iprcparm(mld_slu_ptr_),info)
if (info /= 0) then
ch_err='mld_slu_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
if (info /= 0) then
ch_err='mld_slu_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
class default
info=135
ch_err=a%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end select
call psb_erractionrestore(err_act)
return

89
mlprec/mld_dslud_bld.f90
Unescape Escape View File

@ -74,7 +74,7 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(inout) :: a
type(psb_d_sparse_mat), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(mld_dbaseprec_type), intent(inout) :: p
integer, intent(out) :: info
@ -93,47 +93,56 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (psb_toupper(a%fida) /= 'CSR') then
select type(aa=>a%a)
type is (psb_d_csr_sparse_mat)
!
! WARN: we need to check for a BLOCK distribution (this is the
! distribution required by SuperLU_DIST)
!
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
call psb_loc_to_glob(1,ifrst,desc_a,info)
call psb_loc_to_glob(nrow,ibcheck,desc_a,info)
ibcheck = ibcheck - ifrst + 1
ibcheck = ibcheck - nrow
call psb_amx(ictxt,ibcheck)
if (ibcheck > 0) then
write(0,*) 'Warning: does not look like a BLOCK distribution'
info=135
ch_err = aa%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
endif
mglob = psb_cd_get_global_rows(desc_a)
nzt = aa%get_nzeros()
npr = np
npc = 1
call psb_loc_to_glob(aa%ja(1:nzt),desc_a,info,iact='I')
!
! Compute the LU factorization
!
call mld_dsludist_fact(mglob,nrow,nzt,ifrst,&
& aa%val,aa%irp,aa%ja,p%iprcparm(mld_slud_ptr_),&
& npr, npc, info)
if (info /= 0) then
ch_err='psb_sludist_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
call psb_glob_to_loc(aa%ja(1:nzt),desc_a,info,iact='I')
class default
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif
!
! WARN: we need to check for a BLOCK distribution (this is the
! distribution required by SuperLU_DIST)
!
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
call psb_loc_to_glob(1,ifrst,desc_a,info)
call psb_loc_to_glob(nrow,ibcheck,desc_a,info)
ibcheck = ibcheck - ifrst + 1
ibcheck = ibcheck - nrow
call psb_amx(ictxt,ibcheck)
if (ibcheck > 0) then
write(0,*) 'Warning: does not look like a BLOCK distribution'
endif
mglob = psb_cd_get_global_rows(desc_a)
nzt = psb_sp_get_nnzeros(a)
npr = np
npc = 1
call psb_loc_to_glob(a%ia1(1:nzt),desc_a,info,iact='I')
!
! Compute the LU factorization
!
call mld_dsludist_fact(mglob,nrow,nzt,ifrst,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slud_ptr_),&
& npr, npc, info)
if (info /= 0) then
ch_err='psb_sludist_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
ch_err = aa%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_glob_to_loc(a%ia1(1:nzt),desc_a,info,iact='I')
end select
call psb_erractionrestore(err_act)
return

63
mlprec/mld_dsp_renum.f90
Unescape Escape View File

@ -89,16 +89,18 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a,blck
type(psb_dspmat_type), intent(out) :: atmp
type(psb_d_sparse_mat), intent(in) :: a,blck
type(psb_d_sparse_mat), intent(out) :: atmp
type(mld_dbaseprec_type), intent(inout) :: p
integer, intent(out) :: info
! Local variables
character(len=20) :: name, ch_err
integer :: nztota, nztotb, nztmp, nnr, i,k
integer :: nztota, nztotb, nztmp, nzt2, nnr, i,k, ma, mb
integer, allocatable :: itmp(:), itmp2(:)
integer :: ictxt,np,me, err_act
type(psb_d_coo_sparse_mat) :: cootmp, cootmp2
type(psb_d_csr_sparse_mat) :: csrtmp
real(psb_dpk_) :: t3,t4
if (psb_get_errstatus().ne.0) return
@ -119,18 +121,19 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
! Convert a into the COO format and extend it up to a%m+blck%m rows
! by adding null rows. The converted extended matrix is stored in atmp.
!
nztota=psb_sp_get_nnzeros(a)
nztotb=psb_sp_get_nnzeros(blck)
call psb_spcnv(a,atmp,info,afmt='coo',dupl=psb_dupl_add_)
call psb_rwextd(a%m+blck%m,atmp,info,blck)
nztota=a%get_nzeros()
nztotb=blck%get_nzeros()
ma = a%get_nrows()
mb = blck%get_nrows()
if (p%iprcparm(mld_sub_ren_)==mld_renum_glb_) then
!
! Remember: we have switched IA1=COLS and IA2=ROWS.
! Now identify the set of distinct local column indices.
!
nnr = p%desc_data%matrix_data(psb_n_row_)
nnr = psb_cd_get_local_rows(p%desc_data)
allocate(p%perm(nnr),p%invperm(nnr),itmp2(nnr),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
@ -163,13 +166,16 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
!
! Convert atmp into the CSR format
!
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
nztmp = psb_sp_get_nnzeros(atmp)
call a%cscnv(atmp,info,type='coo',dupl=psb_dupl_add_)
call psb_rwextd(ma+mb,atmp,info,blck)
call atmp%mv_to(csrtmp)
nztmp = csrtmp%get_nzeros()
!
! Realloc the permutation arrays
!
call psb_realloc(atmp%m,p%perm,info)
call psb_realloc(csrtmp%get_nrows(),p%perm,info)
if(info/=0) then
info=4010
ch_err='psb_realloc'
@ -177,7 +183,7 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
goto 9999
end if
call psb_realloc(atmp%m,p%invperm,info)
call psb_realloc(csrtmp%get_nrows(),p%invperm,info)
if(info/=0) then
info=4010
ch_err='psb_realloc'
@ -185,7 +191,7 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
goto 9999
end if
allocate(itmp(max(8,atmp%m+2,nztmp+2)),stat=info)
allocate(itmp(max(8,csrtmp%get_nrows()+2,nztmp+2)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
@ -196,7 +202,7 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
!
! Renumber rows and columns according to the GPS algorithm
!
call gps_reduction(atmp%m,atmp%ia2,atmp%ia1,p%perm,p%invperm,info)
call gps_reduction(csrtmp%get_nrows(),csrtmp%irp,csrtmp%ja,p%perm,p%invperm,info)
if(info/=0) then
info=4010
ch_err='gps_reduction'
@ -207,24 +213,37 @@ subroutine mld_dsp_renum(a,blck,p,atmp,info)
!
! Compute the inverse permutation
!
do k=1, atmp%m
do k=1, csrtmp%get_nrows()
p%invperm(p%perm(k)) = k
enddo
t3 = psb_wtime()
call psb_spcnv(atmp,info,afmt='coo',dupl=psb_dupl_add_)
end if
!
! Rebuild atmp with the new numbering (COO format)
!
nztmp=psb_sp_get_nnzeros(atmp)
call a%cp_to(cootmp)
nztmp=cootmp%get_nzeros()
do i=1,nztmp
atmp%ia1(i) = p%perm(a%ia1(i))
atmp%ia2(i) = p%invperm(a%ia2(i))
cootmp%ia(i) = p%perm(cootmp%ia(i))
cootmp%ja(i) = p%invperm(cootmp%ja(i))
end do
call blck%cp_to(cootmp2)
nzt2 = cootmp2%get_nzeros()
call psb_ensure_size(nztmp+nzt2,cootmp%ia,info)
call psb_ensure_size(nztmp+nzt2,cootmp%ja,info)
call psb_ensure_size(nztmp+nzt2,cootmp%val,info)
do i=1,nzt2
cootmp%ia(nztmp+i) = p%perm(cootmp2%ia(i))
cootmp%ja(nztmp+i) = p%invperm(cootmp2%ja(i))
cootmp%val(nztmp+i) = (cootmp2%val(i))
end do
call psb_spcnv(atmp,info,afmt='coo',dupl=psb_dupl_add_)
call cootmp2%free()
call cootmp%set_nzeros(nztmp+nzt2)
call cootmp%set_dupl(psb_dupl_add_)
call cootmp%fix(info)
call atmp%mv_from(cootmp)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_fixcoo')
goto 9999

8
mlprec/mld_dsub_solve.f90
Unescape Escape View File

@ -200,15 +200,15 @@ subroutine mld_dsub_solve(alpha,prec,x,beta,y,desc_data,trans,work,info)
case('N')
call psb_spsm(done,prec%av(mld_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux)
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if (info == 0) call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,unit='U',choice=psb_none_, work=aux)
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T','C')
call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux)
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if (info == 0) call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,unit='U',choice=psb_none_,work=aux)
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(4001,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999

40
mlprec/mld_dumf_bld.f90
Unescape Escape View File

@ -83,15 +83,14 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
implicit none
! Arguments
type(psb_dspmat_type), intent(inout) :: a
type(psb_d_sparse_mat), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(mld_dbaseprec_type), intent(inout) :: p
integer, intent(out) :: info
! Local variables
integer :: nzt,ictxt,me,np,err_act
integer :: i_err(5)
character(len=20) :: name
integer :: ictxt,me,np,err_act
character(len=20) :: name, ch_err
info=0
name='mld_dumf_bld'
@ -99,27 +98,26 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
if (psb_toupper(a%fida) /= 'CSC') then
info=135
call psb_errpush(info,name,a_err=a%fida)
goto 9999
endif
nzt = psb_sp_get_nnzeros(a)
!
! Compute the LU factorization
!
call mld_dumf_fact(a%m,nzt,&
& a%aspk,a%ia1,a%ia2,&
& p%iprcparm(mld_umf_symptr_),p%iprcparm(mld_umf_numptr_),info)
if (info /= 0) then
i_err(1) = info
info=4110
call psb_errpush(info,name,a_err='mld_umf_fact',i_err=i_err)
select type(aa=>a%a)
!!$ type is (psb_d_csc_sparse_mat)
!!$ call mld_dumf_fact(aa%m,aa%get_nzeros(),&
!!$ & aa%val,aa%ia,aa%icp,&
!!$ & p%iprcparm(mld_umf_symptr_),p%iprcparm(mld_umf_numptr_),info)
!!$
!!$ if (info /= 0) then
!!$ info=4110
!!$ call psb_errpush(info,name,a_err='mld_umf_fact',i_err=(/info,0,0,0,0/))
!!$ goto 9999
!!$ end if
class default
info=135
ch_err = aa%get_fmt()
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end select
call psb_erractionrestore(err_act)
return

450
mlprec/mld_inner_mod.f90
View File

File diff suppressed because it is too large Load Diff

613
mlprec/mld_prec_mod.f90
Unescape Escape View File

@ -48,77 +48,78 @@ module mld_prec_mod
use mld_prec_type
interface mld_precinit
subroutine mld_sprecinit(p,ptype,info,nlev)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
end subroutine mld_sprecinit
!!$ subroutine mld_sprecinit(p,ptype,info,nlev)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ character(len=*), intent(in) :: ptype
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: nlev
!!$ end subroutine mld_sprecinit
subroutine mld_dprecinit(p,ptype,info,nlev)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
end subroutine mld_dprecinit
subroutine mld_cprecinit(p,ptype,info,nlev)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
end subroutine mld_cprecinit
subroutine mld_zprecinit(p,ptype,info,nlev)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: nlev
end subroutine mld_zprecinit
!!$ subroutine mld_cprecinit(p,ptype,info,nlev)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ character(len=*), intent(in) :: ptype
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: nlev
!!$ end subroutine mld_cprecinit
!!$ subroutine mld_zprecinit(p,ptype,info,nlev)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ character(len=*), intent(in) :: ptype
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: nlev
!!$ end subroutine mld_zprecinit
end interface
interface mld_precset
module procedure mld_i_sprecseti, mld_i_sprecsetc, mld_i_sprecsetr,&
& mld_i_dprecseti, mld_i_dprecsetc, mld_i_dprecsetr,&
& mld_i_cprecseti, mld_i_cprecsetc, mld_i_cprecsetr,&
& mld_i_zprecseti, mld_i_zprecsetc, mld_i_zprecsetr
module procedure mld_i_dprecseti, mld_i_dprecsetc, mld_i_dprecsetr
!!$ module procedure mld_i_sprecseti, mld_i_sprecsetc, mld_i_sprecsetr,&
!!$ & mld_i_dprecseti, mld_i_dprecsetc, mld_i_dprecsetr,&
!!$ & mld_i_cprecseti, mld_i_cprecsetc, mld_i_cprecsetr,&
!!$ & mld_i_zprecseti, mld_i_zprecsetc, mld_i_zprecsetr
end interface
interface mld_inner_precset
subroutine mld_sprecseti(p,what,val,info,ilev)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_sprecseti
subroutine mld_sprecsetr(p,what,val,info,ilev)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_sprecsetr
subroutine mld_sprecsetc(p,what,string,info,ilev)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: string
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_sprecsetc
!!$ subroutine mld_sprecseti(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_sprecseti
!!$ subroutine mld_sprecsetr(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_spk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_sprecsetr
!!$ subroutine mld_sprecsetc(p,what,string,info,ilev)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: string
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_sprecsetc
subroutine mld_dprecseti(p,what,val,info,ilev)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -127,7 +128,7 @@ module mld_prec_mod
integer, optional, intent(in) :: ilev
end subroutine mld_dprecseti
subroutine mld_dprecsetr(p,what,val,info,ilev)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -136,7 +137,7 @@ module mld_prec_mod
integer, optional, intent(in) :: ilev
end subroutine mld_dprecsetr
subroutine mld_dprecsetc(p,what,string,info,ilev)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -144,181 +145,181 @@ module mld_prec_mod
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_dprecsetc
subroutine mld_cprecseti(p,what,val,info,ilev)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_cprecseti
subroutine mld_cprecsetr(p,what,val,info,ilev)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_cprecsetr
subroutine mld_cprecsetc(p,what,string,info,ilev)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: string
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_cprecsetc
subroutine mld_zprecseti(p,what,val,info,ilev)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_zprecseti
subroutine mld_zprecsetr(p,what,val,info,ilev)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_zprecsetr
subroutine mld_zprecsetc(p,what,string,info,ilev)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: string
integer, intent(out) :: info
integer, optional, intent(in) :: ilev
end subroutine mld_zprecsetc
!!$ subroutine mld_cprecseti(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_cprecseti
!!$ subroutine mld_cprecsetr(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_spk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_cprecsetr
!!$ subroutine mld_cprecsetc(p,what,string,info,ilev)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: string
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_cprecsetc
!!$ subroutine mld_zprecseti(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_zprecseti
!!$ subroutine mld_zprecsetr(p,what,val,info,ilev)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_dpk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_zprecsetr
!!$ subroutine mld_zprecsetc(p,what,string,info,ilev)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: string
!!$ integer, intent(out) :: info
!!$ integer, optional, intent(in) :: ilev
!!$ end subroutine mld_zprecsetc
end interface
interface mld_precaply
subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
end subroutine mld_sprecaply
subroutine mld_sprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_sprecaply1
subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
end subroutine mld_dprecaply
subroutine mld_dprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_dprecaply1
subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
end subroutine mld_cprecaply
subroutine mld_cprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_cprecaply1
subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
end subroutine mld_zprecaply
subroutine mld_zprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_zprecaply1
end interface
!!$
!!$ interface mld_precaply
!!$ subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_sprec_type), intent(in) :: prec
!!$ real(psb_spk_),intent(in) :: x(:)
!!$ real(psb_spk_),intent(inout) :: y(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ real(psb_spk_),intent(inout), optional, target :: work(:)
!!$ end subroutine mld_sprecaply
!!$ subroutine mld_sprecaply1(prec,x,desc_data,info,trans)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_sprec_type), intent(in) :: prec
!!$ real(psb_spk_),intent(inout) :: x(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ end subroutine mld_sprecaply1
!!$ subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
!!$ use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_dprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_dprec_type), intent(in) :: prec
!!$ real(psb_dpk_),intent(in) :: x(:)
!!$ real(psb_dpk_),intent(inout) :: y(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ real(psb_dpk_),intent(inout), optional, target :: work(:)
!!$ end subroutine mld_dprecaply
!!$ subroutine mld_dprecaply1(prec,x,desc_data,info,trans)
!!$ use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_dprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_dprec_type), intent(in) :: prec
!!$ real(psb_dpk_),intent(inout) :: x(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ end subroutine mld_dprecaply1
!!$ subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_cprec_type), intent(in) :: prec
!!$ complex(psb_spk_),intent(in) :: x(:)
!!$ complex(psb_spk_),intent(inout) :: y(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ complex(psb_spk_),intent(inout), optional, target :: work(:)
!!$ end subroutine mld_cprecaply
!!$ subroutine mld_cprecaply1(prec,x,desc_data,info,trans)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_cprec_type), intent(in) :: prec
!!$ complex(psb_spk_),intent(inout) :: x(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ end subroutine mld_cprecaply1
!!$ subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_zprec_type), intent(in) :: prec
!!$ complex(psb_dpk_),intent(in) :: x(:)
!!$ complex(psb_dpk_),intent(inout) :: y(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ complex(psb_dpk_),intent(inout), optional, target :: work(:)
!!$ end subroutine mld_zprecaply
!!$ subroutine mld_zprecaply1(prec,x,desc_data,info,trans)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(psb_desc_type),intent(in) :: desc_data
!!$ type(mld_zprec_type), intent(in) :: prec
!!$ complex(psb_dpk_),intent(inout) :: x(:)
!!$ integer, intent(out) :: info
!!$ character(len=1), optional :: trans
!!$ end subroutine mld_zprecaply1
!!$ end interface
!!$
interface mld_precbld
subroutine mld_sprecbld(a,desc_a,prec,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
implicit none
type(psb_sspmat_type), intent(in), target :: a
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_sprec_type), intent(inout), target :: prec
integer, intent(out) :: info
!!$ character, intent(in),optional :: upd
end subroutine mld_sprecbld
subroutine mld_dprecbld(a,desc_a,prec,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
implicit none
type(psb_dspmat_type), intent(in), target :: a
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type), intent(inout), target :: prec
integer, intent(out) :: info
!!$ character, intent(in),optional :: upd
end subroutine mld_dprecbld
subroutine mld_cprecbld(a,desc_a,prec,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
implicit none
type(psb_cspmat_type), intent(in), target :: a
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_cprec_type), intent(inout), target :: prec
integer, intent(out) :: info
!!$ character, intent(in),optional :: upd
end subroutine mld_cprecbld
subroutine mld_zprecbld(a,desc_a,prec,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
implicit none
type(psb_zspmat_type), intent(in), target :: a
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_zprec_type), intent(inout) :: prec
integer, intent(out) :: info
@ -328,41 +329,41 @@ module mld_prec_mod
contains
subroutine mld_i_sprecseti(p,what,val,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_sprecseti
subroutine mld_i_sprecsetr(p,what,val,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_sprecsetr
subroutine mld_i_sprecsetc(p,what,val,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_sprec_type
type(mld_sprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_sprecsetc
!!$ subroutine mld_i_sprecseti(p,what,val,info)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_sprecseti
!!$
!!$ subroutine mld_i_sprecsetr(p,what,val,info)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_spk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_sprecsetr
!!$
!!$ subroutine mld_i_sprecsetc(p,what,val,info)
!!$ use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_sprec_type
!!$ type(mld_sprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_sprecsetc
subroutine mld_i_dprecseti(p,what,val,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -373,7 +374,7 @@ contains
end subroutine mld_i_dprecseti
subroutine mld_i_dprecsetr(p,what,val,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -384,7 +385,7 @@ contains
end subroutine mld_i_dprecsetr
subroutine mld_i_dprecsetc(p,what,val,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_dprec_type
type(mld_dprec_type), intent(inout) :: p
integer, intent(in) :: what
@ -394,72 +395,72 @@ contains
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_dprecsetc
subroutine mld_i_cprecseti(p,what,val,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_cprecseti
subroutine mld_i_cprecsetr(p,what,val,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_cprecsetr
subroutine mld_i_cprecsetc(p,what,val,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_prec_type, only : mld_cprec_type
type(mld_cprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_cprecsetc
subroutine mld_i_zprecseti(p,what,val,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_zprecseti
subroutine mld_i_zprecsetr(p,what,val,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_zprecsetr
subroutine mld_i_zprecsetc(p,what,val,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_prec_type, only : mld_zprec_type
type(mld_zprec_type), intent(inout) :: p
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
call mld_inner_precset(p,what,val,info)
end subroutine mld_i_zprecsetc
!!$ subroutine mld_i_cprecseti(p,what,val,info)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_cprecseti
!!$
!!$ subroutine mld_i_cprecsetr(p,what,val,info)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_spk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_cprecsetr
!!$
!!$ subroutine mld_i_cprecsetc(p,what,val,info)
!!$ use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
!!$ use mld_prec_type, only : mld_cprec_type
!!$ type(mld_cprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_cprecsetc
!!$
!!$ subroutine mld_i_zprecseti(p,what,val,info)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ integer, intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_zprecseti
!!$
!!$ subroutine mld_i_zprecsetr(p,what,val,info)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ real(psb_dpk_), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_zprecsetr
!!$
!!$ subroutine mld_i_zprecsetc(p,what,val,info)
!!$ use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
!!$ use mld_prec_type, only : mld_zprec_type
!!$ type(mld_zprec_type), intent(inout) :: p
!!$ integer, intent(in) :: what
!!$ character(len=*), intent(in) :: val
!!$ integer, intent(out) :: info
!!$
!!$ call mld_inner_precset(p,what,val,info)
!!$ end subroutine mld_i_zprecsetc
!!$
end module mld_prec_mod

471
mlprec/mld_prec_type.f90
Unescape Escape View File

@ -64,16 +64,18 @@ module mld_prec_type
! blows up on some systems.
!
use psb_base_mod, only :&
& psb_dspmat_type, psb_zspmat_type,&
& psb_sspmat_type, psb_cspmat_type,&
& psb_d_sparse_mat, psb_z_sparse_mat,&
& psb_s_sparse_mat, psb_c_sparse_mat,&
& psb_desc_type,&
& psb_slinmap_type, psb_dlinmap_type,&
& psb_clinmap_type, psb_zlinmap_type, &
& psb_dpk_, psb_spk_, psb_long_int_k_, &
& psb_sp_free, psb_cdfree, psb_halo_, psb_none_, psb_sum_, psb_avg_, &
& psb_spfree, psb_cdfree, psb_halo_, psb_none_, psb_sum_, psb_avg_, &
& psb_nohalo_, psb_square_root_, psb_toupper, psb_root_,&
& psb_sizeof_int, psb_sizeof_long_int, psb_sizeof_sp, psb_sizeof_dp, psb_sizeof,&
& psb_cd_get_context, psb_info
use psb_prec_mod, only: psb_sprec_type, psb_dprec_type,&
& psb_cprec_type, psb_zprec_type
!
! Type: mld_Tprec_type.
@ -103,9 +105,9 @@ module mld_prec_type
! type(mld_Tbaseprec_type) :: prec
! integer, allocatable :: iprcparm(:)
! real(psb_Tpk_), allocatable :: rprcparm(:)
! type(psb_Tspmat_type) :: ac
! type(psb_T_sparse_mat) :: ac
! type(psb_desc_type) :: desc_ac
! type(psb_Tspmat_type), pointer :: base_a => null()
! type(psb_T_sparse_mat), pointer :: base_a => null()
! type(psb_desc_type), pointer :: base_desc => null()
! type(psb_Tlinmap_type) :: map
! end type mld_Tonelev_type
@ -124,7 +126,7 @@ module mld_prec_type
! desc_ac - type(psb_desc_type).
! The communication descriptor associated to the matrix
! stored in ac.
! base_a - type(psb_zspmat_type), pointer.
! base_a - type(psb_z_sparse_mat), pointer.
! Pointer (really a pointer!) to the local part of the current
! matrix (so we have a unified treatment of residuals).
! We need this to avoid passing explicitly the current matrix
@ -142,7 +144,7 @@ module mld_prec_type
! It holds the smoother (base preconditioner) at a single level.
!
! type mld_Tbaseprec_type
! type(psb_Tspmat_type), allocatable :: av(:)
! type(psb_T_sparse_mat), allocatable :: av(:)
! IntrType(psb_Tpk_), allocatable :: d(:)
! type(psb_desc_type) :: desc_data
! integer, allocatable :: iprcparm(:)
@ -154,7 +156,7 @@ module mld_prec_type
! the kind of the real or complex type, according to the real/complex, single/double
! precision version of MLD2P4.
!
! av - type(psb_Tspmat_type), dimension(:), allocatable(:).
! av - type(psb_T_sparse_mat), dimension(:), allocatable(:).
! The sparse matrices needed to apply the preconditioner at
! the current level ilev.
! av(mld_l_pr_) - The L factor of the ILU factorization of the local
@ -191,7 +193,7 @@ module mld_prec_type
!
type mld_sbaseprec_type
type(psb_sspmat_type), allocatable :: av(:)
type(psb_s_sparse_mat), allocatable :: av(:)
real(psb_spk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
@ -203,19 +205,22 @@ module mld_prec_type
type(mld_sbaseprec_type) :: prec
integer, allocatable :: iprcparm(:)
real(psb_spk_), allocatable :: rprcparm(:)
type(psb_sspmat_type) :: ac
type(psb_s_sparse_mat) :: ac
type(psb_desc_type) :: desc_ac
type(psb_sspmat_type), pointer :: base_a => null()
type(psb_s_sparse_mat), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
type(psb_slinmap_type) :: map
end type mld_sonelev_type
type mld_sprec_type
type, extends(psb_sprec_type) :: mld_sprec_type
type(mld_sonelev_type), allocatable :: precv(:)
contains
procedure, pass(prec) :: s_apply2v => mld_s_apply2v
procedure, pass(prec) :: s_apply1v => mld_s_apply1v
end type mld_sprec_type
type mld_dbaseprec_type
type(psb_dspmat_type), allocatable :: av(:)
type(psb_d_sparse_mat), allocatable :: av(:)
real(psb_dpk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
@ -227,20 +232,23 @@ module mld_prec_type
type(mld_dbaseprec_type) :: prec
integer, allocatable :: iprcparm(:)
real(psb_dpk_), allocatable :: rprcparm(:)
type(psb_dspmat_type) :: ac
type(psb_d_sparse_mat) :: ac
type(psb_desc_type) :: desc_ac
type(psb_dspmat_type), pointer :: base_a => null()
type(psb_d_sparse_mat), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
type(psb_dlinmap_type) :: map
end type mld_donelev_type
type mld_dprec_type
type, extends(psb_dprec_type) :: mld_dprec_type
type(mld_donelev_type), allocatable :: precv(:)
contains
procedure, pass(prec) :: d_apply2v => mld_d_apply2v
procedure, pass(prec) :: d_apply1v => mld_d_apply1v
end type mld_dprec_type
type mld_cbaseprec_type
type(psb_cspmat_type), allocatable :: av(:)
type(psb_c_sparse_mat), allocatable :: av(:)
complex(psb_spk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
@ -252,19 +260,22 @@ module mld_prec_type
type(mld_cbaseprec_type) :: prec
integer, allocatable :: iprcparm(:)
real(psb_spk_), allocatable :: rprcparm(:)
type(psb_cspmat_type) :: ac
type(psb_c_sparse_mat) :: ac
type(psb_desc_type) :: desc_ac
type(psb_cspmat_type), pointer :: base_a => null()
type(psb_c_sparse_mat), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
type(psb_clinmap_type) :: map
end type mld_conelev_type
type mld_cprec_type
type, extends(psb_cprec_type) :: mld_cprec_type
type(mld_conelev_type), allocatable :: precv(:)
contains
procedure, pass(prec) :: c_apply2v => mld_c_apply2v
procedure, pass(prec) :: c_apply1v => mld_c_apply1v
end type mld_cprec_type
type mld_zbaseprec_type
type(psb_zspmat_type), allocatable :: av(:)
type(psb_z_sparse_mat), allocatable :: av(:)
complex(psb_dpk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
@ -276,15 +287,18 @@ module mld_prec_type
type(mld_zbaseprec_type) :: prec
integer, allocatable :: iprcparm(:)
real(psb_dpk_), allocatable :: rprcparm(:)
type(psb_zspmat_type) :: ac
type(psb_z_sparse_mat) :: ac
type(psb_desc_type) :: desc_ac
type(psb_zspmat_type), pointer :: base_a => null()
type(psb_z_sparse_mat), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
type(psb_zlinmap_type) :: map
end type mld_zonelev_type
type mld_zprec_type
type, extends(psb_zprec_type) :: mld_zprec_type
type(mld_zonelev_type), allocatable :: precv(:)
contains
procedure, pass(prec) :: z_apply2v => mld_z_apply2v
procedure, pass(prec) :: z_apply1v => mld_z_apply1v
end type mld_zprec_type
@ -482,6 +496,97 @@ module mld_prec_type
& mld_c_onelev_prec_sizeof, mld_z_onelev_prec_sizeof
end interface
interface mld_precaply
subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
import mld_sprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
end subroutine mld_sprecaply
subroutine mld_sprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_s_sparse_mat, psb_desc_type, psb_spk_
import mld_sprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_sprecaply1
subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
import mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
end subroutine mld_dprecaply
subroutine mld_dprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_d_sparse_mat, psb_desc_type, psb_dpk_
import mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_dprecaply1
subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
import mld_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
end subroutine mld_cprecaply
subroutine mld_cprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_c_sparse_mat, psb_desc_type, psb_spk_
import mld_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_cprecaply1
subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
import mld_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
end subroutine mld_zprecaply
subroutine mld_zprecaply1(prec,x,desc_data,info,trans)
use psb_base_mod, only : psb_z_sparse_mat, psb_desc_type, psb_dpk_
import mld_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine mld_zprecaply1
end interface
contains
!
@ -1886,12 +1991,7 @@ contains
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
! return
end if
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
end if
@ -1937,7 +2037,7 @@ contains
! for the inner UMFPACK or SLU stuff
call mld_precfree(p%prec,info)
call psb_sp_free(p%ac,info)
call p%ac%free()
if (allocated(p%desc_ac%matrix_data)) &
& call psb_cdfree(p%desc_ac,info)
@ -1996,12 +2096,7 @@ contains
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
! return
end if
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
end if
@ -2052,7 +2147,7 @@ contains
! for the inner UMFPACK or SLU stuff
call mld_precfree(p%prec,info)
call psb_sp_free(p%ac,info)
call p%ac%free()
if (allocated(p%desc_ac%matrix_data)) &
& call psb_cdfree(p%desc_ac,info)
@ -2105,12 +2200,7 @@ contains
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
! return
end if
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
@ -2154,7 +2244,7 @@ contains
! for the inner UMFPACK or SLU stuff
call mld_precfree(p%prec,info)
call psb_sp_free(p%ac,info)
call p%ac%free()
if (allocated(p%desc_ac%matrix_data)) &
& call psb_cdfree(p%desc_ac,info)
@ -2207,12 +2297,7 @@ contains
if (allocated(p%av)) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
! Actually, we don't care here about this.
! Just let it go.
! return
end if
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
@ -2260,7 +2345,7 @@ contains
! for the inner UMFPACK or SLU stuff
call mld_precfree(p%prec,info)
call psb_sp_free(p%ac,info)
call p%ac%free()
if (allocated(p%desc_ac%matrix_data)) &
& call psb_cdfree(p%desc_ac,info)
@ -2482,4 +2567,286 @@ contains
end subroutine mld_zprec_free
subroutine mld_s_apply2v(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='s_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_sprec_type)
!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_s_apply2v
subroutine mld_s_apply1v(prec,x,desc_data,info,trans)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_sprec_type), intent(in) :: prec
real(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
Integer :: err_act
character(len=20) :: name='s_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_sprec_type)
!!$ call mld_precaply(prec,x,desc_data,info,trans)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_s_apply1v
subroutine mld_d_apply2v(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='d_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
type is (mld_dprec_type)
call mld_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_d_apply2v
subroutine mld_d_apply1v(prec,x,desc_data,info,trans)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
Integer :: err_act
character(len=20) :: name='d_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
type is (mld_dprec_type)
call mld_precaply(prec,x,desc_data,info,trans)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_d_apply1v
subroutine mld_c_apply2v(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='s_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_cprec_type)
!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_c_apply2v
subroutine mld_c_apply1v(prec,x,desc_data,info,trans)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
Integer :: err_act
character(len=20) :: name='c_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_cprec_type)
!!$ call mld_precaply(prec,x,desc_data,info,trans)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_c_apply1v
subroutine mld_z_apply2v(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='z_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_zprec_type)
!!$ call mld_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_z_apply2v
subroutine mld_z_apply1v(prec,x,desc_data,info,trans)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
Integer :: err_act
character(len=20) :: name='z_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
!!$ type is (mld_zprec_type)
!!$ call mld_precaply(prec,x,desc_data,info,trans)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_z_apply1v
end module mld_prec_type

4
tests/pdegen/Makefile
Unescape Escape View File

@ -2,8 +2,8 @@ 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
MLD_LIB=-L$(MLDLIBDIR) -lmld_prec
PSBLAS_LIB= -L$(PSBDIR) -lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base
FINCLUDES=$(FMFLAG). $(FMFLAG)$(MLDLIBDIR) $(FMFLAG)$(PSBDIR) $(FIFLAG).

9
tests/pdegen/ppde.f90
Unescape Escape View File

@ -86,8 +86,8 @@ program ppde
real(psb_dpk_) :: t1, t2, tprec
! sparse matrix and preconditioner
type(psb_dspmat_type) :: a
type(mld_dprec_type) :: prec
type(psb_d_sparse_mat) :: a
type(mld_dprec_type) :: prec
! descriptor
type(psb_desc_type) :: desc_a
! dense matrices
@ -417,7 +417,7 @@ contains
type(psb_desc_type) :: desc_a
integer :: ictxt, info
character :: afmt*5
type(psb_dspmat_type) :: a
type(psb_d_sparse_mat) :: a
real(psb_dpk_) :: zt(nb),glob_x,glob_y,glob_z
integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k
integer :: x,y,z,ia,indx_owner
@ -679,8 +679,9 @@ contains
call psb_amx(ictxt,tasb)
call psb_amx(ictxt,ttot)
if(iam == psb_root_) then
ch_err = a%get_fmt()
write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
& a%fida(1:3)
& ch_err(1:3)
write(*,'("-allocation time : ",es12.5)') talc
write(*,'("-coeff. gen. time : ",es12.5)') tgen
write(*,'("-assembly time : ",es12.5)') tasb

Write Preview
Loading…
Cancel
Save