Merge branch 'master' into maint-1.0
Salvatore Filippone
4 years ago
commit
636600f1c7
@ -0,0 +1,120 @@
|
|||||||
|
##########################################################
|
||||||
|
.mod=@MODEXT@
|
||||||
|
.fh=.fh
|
||||||
|
.SUFFIXES:
|
||||||
|
.SUFFIXES: .f90 .F90 .f .F .c .cpp .o
|
||||||
|
# The following ones are the variables used by the PSBLAS make scripts.
|
||||||
|
|
||||||
|
FC=@FC@
|
||||||
|
CC=@CC@
|
||||||
|
CXX=@CXX@
|
||||||
|
FCOPT=@FCOPT@
|
||||||
|
CCOPT=@CCOPT@
|
||||||
|
CXXOPT=@CXXOPT@
|
||||||
|
FMFLAG=@FMFLAG@
|
||||||
|
FIFLAG=@FIFLAG@
|
||||||
|
EXTRA_OPT=@EXTRA_OPT@
|
||||||
|
|
||||||
|
# These three should be always set!
|
||||||
|
MPFC=@MPIFC@
|
||||||
|
MPCC=@MPICC@
|
||||||
|
MPCXX=@MPICXX@
|
||||||
|
|
||||||
|
FLINK=@FLINK@
|
||||||
|
|
||||||
|
LIBS=@LIBS@
|
||||||
|
|
||||||
|
# BLAS, BLACS and METIS libraries.
|
||||||
|
BLAS=@BLAS_LIBS@
|
||||||
|
METIS_LIB=@METIS_LIBS@
|
||||||
|
LAPACK=@LAPACK_LIBS@
|
||||||
|
|
||||||
|
PSBFDEFINES=@FDEFINES@
|
||||||
|
PSBCDEFINES=@CDEFINES@
|
||||||
|
PSBCXXDEFINES=@CDEFINES@
|
||||||
|
AR=@AR@
|
||||||
|
RANLIB=@RANLIB@
|
||||||
|
|
||||||
|
##########################################################
|
||||||
|
# #
|
||||||
|
# Note: directories external to the AMG4PSBLAS subtree #
|
||||||
|
# must be specified here with absolute pathnames #
|
||||||
|
# #
|
||||||
|
##########################################################
|
||||||
|
PSBLASDIR=@PSBLAS_DIR@
|
||||||
|
PSBLAS_INCDIR=@PSBLAS_INCDIR@
|
||||||
|
PSBLAS_MODDIR=@PSBLAS_MODDIR@
|
||||||
|
PSBLAS_LIBDIR=@PSBLAS_LIBDIR@
|
||||||
|
PSBLAS_INCLUDES=@PSBLAS_INCLUDES@
|
||||||
|
PSBLAS_LIBS=@PSBLAS_LIBS@
|
||||||
|
PSBBASEMODNAME=psb_base_mod
|
||||||
|
PSBPRECMODNAME=psb_prec_mod
|
||||||
|
PSBMETHDMODNAME=psb_krylov_mod
|
||||||
|
PSBUTILMODNAME=psb_util_mod
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
INSTALL=@INSTALL@
|
||||||
|
INSTALL_DATA=@INSTALL_DATA@
|
||||||
|
INSTALL_DIR=@INSTALL_DIR@
|
||||||
|
INSTALL_LIBDIR=@INSTALL_LIBDIR@
|
||||||
|
INSTALL_INCLUDEDIR=@INSTALL_INCLUDEDIR@
|
||||||
|
INSTALL_MODULESDIR=@INSTALL_MODULESDIR@
|
||||||
|
INSTALL_DOCSDIR=@INSTALL_DOCSDIR@
|
||||||
|
INSTALL_SAMPLESDIR=@INSTALL_SAMPLESDIR@
|
||||||
|
|
||||||
|
|
||||||
|
##########################################################
|
||||||
|
# #
|
||||||
|
# Additional defines and libraries for multilevel #
|
||||||
|
# Note that these libraries should be compatible #
|
||||||
|
# (compiled with) the compilers specified in the #
|
||||||
|
# PSBLAS main Make.inc #
|
||||||
|
# #
|
||||||
|
# Examples: #
|
||||||
|
# MUMPSLIBS=-ldmumps -lmumps_common #
|
||||||
|
# -lpord -L/path/to/MUMPS/lib #
|
||||||
|
# MUMPSFLAGS=-DHave_MUMPS_ -I/path/to/MUMPS/include #
|
||||||
|
# #
|
||||||
|
# UMFLIBS=-lumfpack -lamd -L/path/to/UMFPACK #
|
||||||
|
# UMFFLAGS=-DHave_UMF_ -I/path/to/UMFPACK #
|
||||||
|
# #
|
||||||
|
# SLULIBS=-lslu -L/path/to/SuperLU #
|
||||||
|
# SLUFLAGS=-DHave_SLU_ -I/path/to/SuperLU #
|
||||||
|
# #
|
||||||
|
# SLUDISTLIBS=-lslud -L/path/to/SuperLUDist #
|
||||||
|
# SLUDISTFLAGS=-DHave_SLUDist_ -I/path/to/SuperLUDist #
|
||||||
|
# #
|
||||||
|
##########################################################
|
||||||
|
|
||||||
|
MUMPSLIBS=@MUMPS_LIBS@
|
||||||
|
MUMPSFLAGS=@MUMPS_FLAGS@
|
||||||
|
|
||||||
|
SLULIBS=@SLU_LIBS@
|
||||||
|
SLUFLAGS=@SLU_FLAGS@
|
||||||
|
|
||||||
|
SLUDISTLIBS=@SLUDIST_LIBS@
|
||||||
|
SLUDISTFLAGS=@SLUDIST_FLAGS@
|
||||||
|
|
||||||
|
UMFLIBS=@UMF_LIBS@
|
||||||
|
UMFFLAGS=@UMF_FLAGS@
|
||||||
|
|
||||||
|
EXTRALIBS=@EXTRA_LIBS@
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
@COMPILERULES@
|
||||||
|
|
||||||
|
#
|
||||||
|
AMGCDEFINES=$(MUMPSFLAGS) $(SLUFLAGS) $(UMFFLAGS) $(SLUDISTFLAGS) $(PSBCDEFINES)
|
||||||
|
CDEFINES=$(AMGCDEFINES)
|
||||||
|
AMGFDEFINES=@AMGFDEFINES@ $(PSBFDEFINES)
|
||||||
|
FDEFINES=$(AMGFDEFINES)
|
||||||
|
|
||||||
|
CXXDEFINES=@AMGCXXDEFINES@ $(PSBCXXDEFINES)
|
||||||
|
|
||||||
|
AMGLDLIBS=$(MUMPSLIBS) $(SLULIBS) $(SLUDISTLIBS) $(UMFLIBS) $(EXTRALIBS) $(PSBLDLIBS) -lstdc++
|
||||||
|
LDLIBS=$(AMGLDLIBS)
|
||||||
|
|
||||||
|
|
||||||
@ -0,0 +1,26 @@
|
|||||||
|
#!/bin/sh
|
||||||
|
# $Id$
|
||||||
|
|
||||||
|
# This file is still not complete.
|
||||||
|
|
||||||
|
# this should create configure from configure.ac..
|
||||||
|
|
||||||
|
#touch NEWS AUTHORS ChangeLog COPYING INSTALL missing install-sh
|
||||||
|
|
||||||
|
#libtoolize
|
||||||
|
|
||||||
|
aclocal -I config || { echo "no aclocal ?" ; exit 1 ; }
|
||||||
|
|
||||||
|
# we produce a configure script
|
||||||
|
autoconf -o configure_n configure_n.ac || { echo "no autoconf ?" ; exit 1 ; }
|
||||||
|
|
||||||
|
# we produce a brand new Makefile
|
||||||
|
#automake -i --add-missing || { echo "no automake ?" ; exit 1 ; }
|
||||||
|
#automake -i || { echo "no automake ?" ; exit 1 ; }
|
||||||
|
|
||||||
|
# The automake required for autogen.sh'in this package is 1.10#.
|
||||||
|
# So users who want to build from the svn repository are required to use this version.
|
||||||
|
#
|
||||||
|
# Users building from the tarball shouldn't bother, of course,
|
||||||
|
# because they get the configure script generated from the tarball maintainer.
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,991 @@
|
|||||||
|
dnl $Id$
|
||||||
|
|
||||||
|
dnl Process this file with autoconf to produce a configure script.
|
||||||
|
dnl
|
||||||
|
dnl usage : aclocal -I config/ && autoconf && ./configure && make
|
||||||
|
|
||||||
|
dnl then : VAR=VAL ./configure
|
||||||
|
|
||||||
|
dnl In some configurations (AIX) the next line is needed:
|
||||||
|
dnl MPIFC=mpxlf95 ./configure
|
||||||
|
|
||||||
|
dnl then : ./configure VAR=VAL
|
||||||
|
dnl then : ./configure --help=short
|
||||||
|
dnl then : ./configure --help
|
||||||
|
|
||||||
|
dnl the PSBLAS modules get this task difficult to accomplish!
|
||||||
|
dnl SEE : --module-path --include-path
|
||||||
|
|
||||||
|
dnl NOTE : There is no cross compilation support.
|
||||||
|
|
||||||
|
dnl NOTE : missing ifort and kl* library handling..
|
||||||
|
dnl NOTE : odd configurations like ifc + gcc still await in the mist of the unknown
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
###############################################################################
|
||||||
|
#
|
||||||
|
# This script is used by the PSBLAS to determine the compilers, linkers, and
|
||||||
|
# libraries to build its libraries executable code.
|
||||||
|
# Its behaviour is driven on the compiler it finds or it is dictated to work
|
||||||
|
# with.
|
||||||
|
#
|
||||||
|
###############################################################################
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
# NOTE: the literal for version (the second argument to AC_INIT should be a literal!)
|
||||||
|
AC_INIT([AMG4PSBLAS],1.0.0, [https://github.com/sfilippone/amg4psblas/issues])
|
||||||
|
|
||||||
|
# VERSION is the file containing the PSBLAS version code
|
||||||
|
# FIXME
|
||||||
|
amg4psblas_cv_version="1.0.0"
|
||||||
|
|
||||||
|
# A sample source file
|
||||||
|
AC_CONFIG_SRCDIR([amgprec/amg_prec_type.f90])
|
||||||
|
# Our custom M4 macros are in the 'config' directory
|
||||||
|
AC_CONFIG_MACRO_DIR([config])
|
||||||
|
AC_MSG_NOTICE([
|
||||||
|
--------------------------------------------------------------------------------
|
||||||
|
Welcome to the $PACKAGE_NAME $amg4psblas_cv_version configure Script.
|
||||||
|
|
||||||
|
This creates Make.inc, but if you read carefully the
|
||||||
|
documentation, you can make your own by hand for your needs.
|
||||||
|
|
||||||
|
./configure --with-psblas=/path/to/psblas
|
||||||
|
See ./configure --help=short fore more info.
|
||||||
|
--------------------------------------------------------------------------------
|
||||||
|
])
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# FLAGS and LIBS user customization
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
dnl NOTE : no spaces before the comma, and no brackets before the second argument!
|
||||||
|
PAC_ARG_WITH_PSBLAS
|
||||||
|
PSBLAS_DIR="$pac_cv_psblas_dir";
|
||||||
|
PSBLAS_INCDIR="$pac_cv_psblas_incdir";
|
||||||
|
PSBLAS_MODDIR="$pac_cv_psblas_moddir";
|
||||||
|
PSBLAS_LIBDIR="$pac_cv_psblas_libdir";
|
||||||
|
AC_MSG_CHECKING([for PSBLAS install dir])
|
||||||
|
if test "X$PSBLAS_DIR" != "X" ; then
|
||||||
|
case $PSBLAS_DIR in
|
||||||
|
/*) ;;
|
||||||
|
*) AC_MSG_ERROR([The PSBLAS installation dir must be an absolute pathname
|
||||||
|
specified with --with-psblas=/path/to/psblas])
|
||||||
|
esac
|
||||||
|
if test ! -d "$PSBLAS_DIR" ; then
|
||||||
|
AC_MSG_ERROR([Could not find PSBLAS build dir $PSBLAS_DIR!])
|
||||||
|
fi
|
||||||
|
AC_MSG_RESULT([$PSBLAS_DIR])
|
||||||
|
fi
|
||||||
|
|
||||||
|
AM_INIT_AUTOMAKE
|
||||||
|
dnl Specify required version of autoconf.
|
||||||
|
AC_PREREQ(2.59)
|
||||||
|
#
|
||||||
|
# Installation.
|
||||||
|
#
|
||||||
|
#
|
||||||
|
AC_PROG_INSTALL
|
||||||
|
|
||||||
|
AC_MSG_CHECKING([where to install])
|
||||||
|
case $prefix in
|
||||||
|
\/* ) eval "INSTALL_DIR=$prefix";;
|
||||||
|
* ) eval "INSTALL_DIR=/usr/local/amg4psblas";;
|
||||||
|
esac
|
||||||
|
case $libdir in
|
||||||
|
\/* ) eval "INSTALL_LIBDIR=$libdir";;
|
||||||
|
* ) eval "INSTALL_LIBDIR=$INSTALL_DIR/lib";;
|
||||||
|
esac
|
||||||
|
case $includedir in
|
||||||
|
\/* ) eval "INSTALL_INCLUDEDIR=$includedir";;
|
||||||
|
* ) eval "INSTALL_INCLUDEDIR=$INSTALL_DIR/include";;
|
||||||
|
esac
|
||||||
|
INSTALL_MODULESDIR=$INSTALL_DIR/modules
|
||||||
|
case $docsdir in
|
||||||
|
\/* ) eval "INSTALL_DOCSDIR=$docsdir";;
|
||||||
|
* ) eval "INSTALL_DOCSDIR=$INSTALL_DIR/docs";;
|
||||||
|
esac
|
||||||
|
case $samplesdir in
|
||||||
|
\/* ) eval "INSTALL_SAMPLESDIR=$samplesdir";;
|
||||||
|
* ) eval "INSTALL_SAMPLESDIR=$INSTALL_DIR/samples";;
|
||||||
|
esac
|
||||||
|
AC_MSG_RESULT([$INSTALL_DIR $INSTALL_INCLUDEDIR $INSTALL_MODULESDIR $INSTALL_LIBDIR $INSTALL_DOCSDIR $INSTALL_SAMPLESDIR])
|
||||||
|
|
||||||
|
dnl
|
||||||
|
dnl We set our own FC flags, ignore those from AC_PROG_FC but not those from the
|
||||||
|
dnl environment variable. Same for C
|
||||||
|
dnl
|
||||||
|
save_FCFLAGS="$FCFLAGS";
|
||||||
|
AC_PROG_FC([ftn xlf2003_r xlf2003 xlf95_r xlf95 xlf90 xlf pgf95 pgf90 ifort ifc nagfor gfortran])
|
||||||
|
FCFLAGS="$save_FCFLAGS";
|
||||||
|
save_CFLAGS="$CFLAGS";
|
||||||
|
AC_PROG_CC([cc xlc pgcc icc gcc ])
|
||||||
|
CFLAGS="$save_CFLAGS";
|
||||||
|
save_CXXFLAGS="$CXXFLAGS";
|
||||||
|
AC_PROG_CXX([CC xlc++ icpc g++])
|
||||||
|
CXXFLAGS="$save_CXXFLAGS";
|
||||||
|
dnl AC_PROG_CXX
|
||||||
|
|
||||||
|
dnl AC_PROG_F90 doesn't exist, at the time of writing this !
|
||||||
|
dnl AC_PROG_F90
|
||||||
|
|
||||||
|
# Sanity checks, although redundant (useful when debugging this configure.ac)!
|
||||||
|
if test "X$FC" == "X" ; then
|
||||||
|
AC_MSG_ERROR([Problem : No Fortran compiler specified nor found!])
|
||||||
|
fi
|
||||||
|
|
||||||
|
if eval "$FC -qversion 2>&1 | grep XL 2>/dev/null" ; then
|
||||||
|
# Some configurations of the XLF want "-WF," prepended to -D.. flags.
|
||||||
|
# TODO : discover the exact conditions when the usage of -WF is needed.
|
||||||
|
amg_cv_define_prepend="-WF,"
|
||||||
|
if eval "$MPIFC -qversion 2>&1 | grep -e\"Version: 10\.\" 2>/dev/null"; then
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DXLF_10 $FDEFINES"
|
||||||
|
fi
|
||||||
|
|
||||||
|
# Note : there could be problems with old xlf compiler versions ( <10.1 )
|
||||||
|
# since (as far as it is known to us) -WF, is not used in earlier versions.
|
||||||
|
# More problems could be undocumented yet.
|
||||||
|
fi
|
||||||
|
|
||||||
|
if test "X$CC" == "X" ; then
|
||||||
|
AC_MSG_ERROR([Problem : No C compiler specified nor found!])
|
||||||
|
fi
|
||||||
|
AC_PROG_CC_STDC()
|
||||||
|
if test "x$ac_cv_prog_cc_stdc" == "xno" ; then
|
||||||
|
AC_MSG_ERROR([Problem : Need a C99 compiler ! ])
|
||||||
|
else
|
||||||
|
C99OPT="$ac_cv_prog_cc_stdc";
|
||||||
|
fi
|
||||||
|
###############################################################################
|
||||||
|
# Suitable MPI compilers detection
|
||||||
|
###############################################################################
|
||||||
|
# Note: Someday we will contemplate a fake MPI - configured version of PSBLAS
|
||||||
|
###############################################################################
|
||||||
|
# First check whether the user required our serial (fake) mpi.
|
||||||
|
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";
|
||||||
|
MPICC="$CC";
|
||||||
|
MPICXX="$CXX";
|
||||||
|
CXXDEFINES="-DSERIAL_MPI $CXXDEFINES";
|
||||||
|
|
||||||
|
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 mpiicc mpcc mpicc cc])
|
||||||
|
fi
|
||||||
|
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for C]])])
|
||||||
|
|
||||||
|
|
||||||
|
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],[mpxlf2003_r mpxlf2003 mpxlf95_r mpxlf90 mpiifort mpf95 mpf90 mpifort mpif95 mpif90 ftn ])
|
||||||
|
fi
|
||||||
|
|
||||||
|
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for Fortran]])])
|
||||||
|
AC_LANG([C++])
|
||||||
|
if test "X$MPICXX" = "X" ; then
|
||||||
|
# This is our MPICC compiler preference: it will override ACX_MPI's first try.
|
||||||
|
AC_CHECK_PROGS([MPICXX],[mpxlc++ mpiicpc mpicxx])
|
||||||
|
fi
|
||||||
|
ACX_MPI([], [AC_MSG_ERROR([[Cannot find any suitable MPI implementation for C++]])])
|
||||||
|
AC_LANG([Fortran])
|
||||||
|
|
||||||
|
FC="$MPIFC" ;
|
||||||
|
CC="$MPICC";
|
||||||
|
CXX="$MPICXX";
|
||||||
|
fi
|
||||||
|
|
||||||
|
AC_LANG([C])
|
||||||
|
|
||||||
|
dnl Now on, MPIFC should be set, and MPICC
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Sanity checks, although redundant (useful when debugging this configure.ac)!
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
if test "X$MPIFC" == "X" ; then
|
||||||
|
AC_MSG_ERROR([Problem : No MPI Fortran compiler specified nor found!])
|
||||||
|
fi
|
||||||
|
|
||||||
|
if test "X$MPICC" == "X" ; then
|
||||||
|
AC_MSG_ERROR([Problem : No MPI C compiler specified nor found!])
|
||||||
|
fi
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# FLAGS and LIBS user customization
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
dnl NOTE : no spaces before the comma, and no brackets before the second argument!
|
||||||
|
PAC_ARG_WITH_FLAGS(ccopt,CCOPT)
|
||||||
|
PAC_ARG_WITH_FLAGS(cxxopt,CXXOPT)
|
||||||
|
PAC_ARG_WITH_FLAGS(fcopt,FCOPT)
|
||||||
|
PAC_ARG_WITH_LIBS
|
||||||
|
PAC_ARG_WITH_FLAGS(clibs,CLIBS)
|
||||||
|
PAC_ARG_WITH_FLAGS(flibs,FLIBS)
|
||||||
|
PAC_ARG_WITH_FLAGS(cxxlibs,CXXLIBS)
|
||||||
|
|
||||||
|
dnl candidates for removal:
|
||||||
|
PAC_ARG_WITH_FLAGS(library-path,LIBRARYPATH)
|
||||||
|
PAC_ARG_WITH_FLAGS(include-path,INCLUDEPATH)
|
||||||
|
PAC_ARG_WITH_FLAGS(module-path,MODULE_PATH)
|
||||||
|
|
||||||
|
# we just gave the user the chance to append values to these variables
|
||||||
|
PAC_ARG_WITH_EXTRA_LIBS
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Sanity checks, although redundant (useful when debugging this configure.ac)!
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Compiler identification (sadly, it is necessary)
|
||||||
|
###############################################################################
|
||||||
|
psblas_cv_fc=""
|
||||||
|
|
||||||
|
dnl Do we use gfortran & co ? Compiler identification.
|
||||||
|
dnl NOTE : in /autoconf/autoconf/fortran.m4 there are plenty of better tests!
|
||||||
|
|
||||||
|
PAC_CHECK_HAVE_GFORTRAN(
|
||||||
|
[psblas_cv_fc="gcc"],
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_CHECK_HAVE_CRAYFTN(
|
||||||
|
[psblas_cv_fc="cray"],
|
||||||
|
)
|
||||||
|
|
||||||
|
if test x"$psblas_cv_fc" == "x" ; then
|
||||||
|
if eval "$MPIFC -qversion 2>&1 | grep XL 2>/dev/null" ; then
|
||||||
|
psblas_cv_fc="xlf"
|
||||||
|
# Some configurations of the XLF want "-WF," prepended to -D.. flags.
|
||||||
|
# TODO : discover the exact conditions when the usage of -WF is needed.
|
||||||
|
psblas_cv_define_prepend="-WF,"
|
||||||
|
if eval "$MPIFC -qversion 2>&1 | grep -e\"Version: 10\.\" 2>/dev/null"; then
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DXLF_10 $FDEFINES"
|
||||||
|
fi
|
||||||
|
|
||||||
|
# Note : there could be problems with old xlf compiler versions ( <10.1 )
|
||||||
|
# since (as far as it is known to us) -WF, is not used in earlier versions.
|
||||||
|
# More problems could be undocumented yet.
|
||||||
|
elif eval "$MPIFC -V 2>&1 | grep Sun 2>/dev/null" ; then
|
||||||
|
# Sun compiler detection
|
||||||
|
|
||||||
|
psblas_cv_fc="sun"
|
||||||
|
elif eval "$MPIFC -V 2>&1 | grep Portland 2>/dev/null" ; then
|
||||||
|
# Portland group compiler detection
|
||||||
|
|
||||||
|
psblas_cv_fc="pg"
|
||||||
|
elif eval "$MPIFC -V 2>&1 | grep Intel.*Fortran.*Compiler 2>/dev/null" ; then
|
||||||
|
# Intel compiler identification
|
||||||
|
|
||||||
|
psblas_cv_fc="ifc"
|
||||||
|
elif eval "$MPIFC -v 2>&1 | grep NAG 2>/dev/null" ; then
|
||||||
|
psblas_cv_fc="nag"
|
||||||
|
FC="$MPIFC"
|
||||||
|
else
|
||||||
|
psblas_cv_fc=""
|
||||||
|
# unsupported MPI Fortran compiler
|
||||||
|
AC_MSG_NOTICE([[Unknown Fortran compiler, proceeding with fingers crossed !]])
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
if test "X$psblas_cv_fc" == "Xgcc" ; then
|
||||||
|
PAC_HAVE_MODERN_GFORTRAN(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Bailing out.])]
|
||||||
|
)
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Linking, symbol mangling, and misc tests
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
# Note : This is functional to Make.inc rules and structure (see below).
|
||||||
|
AC_LANG([C])
|
||||||
|
AC_CHECK_SIZEOF(void *)
|
||||||
|
# Define for platforms with 64 bit (void * ) pointers
|
||||||
|
if test X"$ac_cv_sizeof_void_p" == X"8" ; then
|
||||||
|
CDEFINES="-DPtr64Bits $CDEFINES"
|
||||||
|
fi
|
||||||
|
AC_LANG([Fortran])
|
||||||
|
__AC_FC_NAME_MANGLING
|
||||||
|
if test "X$psblas_cv_fc" == X"pg" ; then
|
||||||
|
FC=$save_FC
|
||||||
|
fi
|
||||||
|
AC_LANG([C])
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling: $ac_cv_fc_mangling])
|
||||||
|
[pac_fc_case=${ac_cv_fc_mangling%%,*}]
|
||||||
|
[pac_fc_under=${ac_cv_fc_mangling#*,}]
|
||||||
|
[pac_fc_sec_under=${pac_fc_under#*,}]
|
||||||
|
[pac_fc_sec_under=${pac_fc_sec_under# }]
|
||||||
|
[pac_fc_under=${pac_fc_under%%,*}]
|
||||||
|
[pac_fc_under=${pac_fc_under# }]
|
||||||
|
AC_MSG_CHECKING([defines for C/Fortran name interfaces])
|
||||||
|
if test "x$pac_fc_case" == "xlower case"; then
|
||||||
|
if test "x$pac_fc_under" == "xunderscore"; then
|
||||||
|
if test "x$pac_fc_sec_under" == "xno extra underscore"; then
|
||||||
|
pac_f_c_names="-DLowerUnderscore"
|
||||||
|
elif test "x$pac_fc_sec_under" == "xextra underscore"; then
|
||||||
|
pac_f_c_names="-DLowerDoubleUnderscore"
|
||||||
|
else
|
||||||
|
pac_f_c_names="-DUNKNOWN"
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling extra underscore unknown case])
|
||||||
|
fi
|
||||||
|
elif test "x$pac_fc_under" == "xno underscore"; then
|
||||||
|
pac_f_c_names="-DLowerCase"
|
||||||
|
else
|
||||||
|
pac_f_c_names="-DUNKNOWN"
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling underscore unknown case])
|
||||||
|
fi
|
||||||
|
elif test "x$pac_fc_case" == "xupper case"; then
|
||||||
|
if test "x$pac_fc_under" == "xunderscore"; then
|
||||||
|
if test "x$pac_fc_sec_under" == "xno extra underscore"; then
|
||||||
|
pac_f_c_names="-DUpperUnderscore"
|
||||||
|
elif test "x$pac_fc_sec_under" == "xextra underscore"; then
|
||||||
|
pac_f_c_names="-DUpperDoubleUnderscore"
|
||||||
|
else
|
||||||
|
pac_f_c_names="-DUNKNOWN"
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling extra underscore unknown case])
|
||||||
|
fi
|
||||||
|
elif test "x$pac_fc_under" == "xno underscore"; then
|
||||||
|
pac_f_c_names="-DUpperCase"
|
||||||
|
else
|
||||||
|
pac_f_c_names="-DUNKNOWN"
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling underscore unknown case])
|
||||||
|
fi
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling UPPERCASE not handled])
|
||||||
|
else
|
||||||
|
pac_f_c_names="-DUNKNOWN"
|
||||||
|
dnl AC_MSG_NOTICE([Fortran name mangling unknown case])
|
||||||
|
fi
|
||||||
|
CDEFINES="$pac_f_c_names $CDEFINES"
|
||||||
|
|
||||||
|
AC_MSG_RESULT([ $pac_f_c_names ])
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Make.inc generation logic
|
||||||
|
###############################################################################
|
||||||
|
# Honor CFLAGS if they were specified explicitly, but --with-ccopt take precedence
|
||||||
|
if test "X$CCOPT" == "X" ; then
|
||||||
|
CCOPT="$CFLAGS";
|
||||||
|
fi
|
||||||
|
if test "X$CCOPT" == "X" ; then
|
||||||
|
if test "X$psblas_cv_fc" == "Xgcc" ; then
|
||||||
|
# note that no space should be placed around the equality symbol in assignements
|
||||||
|
# Note : 'native' is valid _only_ on GCC/x86 (32/64 bits)
|
||||||
|
CCOPT="-O3 $CCOPT"
|
||||||
|
|
||||||
|
elif test "X$psblas_cv_fc" == X"xlf" ; then
|
||||||
|
# XL compiler : consider using -qarch=auto
|
||||||
|
CCOPT="-O3 -qarch=auto $CCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"ifc" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CCOPT="-O3 $CCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"pg" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CCOPT="-fast $CCOPT"
|
||||||
|
# NOTE : PG & Sun use -fast instead -O3
|
||||||
|
elif test "X$psblas_cv_fc" == X"sun" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CCOPT="-fast $CCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"cray" ; then
|
||||||
|
CCOPT="-O3 $CCOPT"
|
||||||
|
MPICC="cc"
|
||||||
|
elif test "X$psblas_cv_fc" == X"nag" ; then
|
||||||
|
# using GCC in conjunction with NAG.
|
||||||
|
CCOPT="-O2"
|
||||||
|
else
|
||||||
|
CCOPT="-O2 $CCOPT"
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
#CFLAGS="${CCOPT}"
|
||||||
|
if test "X$CXXOPT" == "X" ; then
|
||||||
|
CXXOPT="$CXXFLAGS";
|
||||||
|
fi
|
||||||
|
if test "X$CXXOPT" == "X" ; then
|
||||||
|
if test "X$psblas_cv_fc" == "Xgcc" ; then
|
||||||
|
# note that no space should be placed around the equality symbol in assignements
|
||||||
|
# Note : 'native' is valid _only_ on GCC/x86 (32/64 bits)
|
||||||
|
CXXOPT="-g -O3 $CXXOPT"
|
||||||
|
|
||||||
|
elif test "X$psblas_cv_fc" == X"xlf" ; then
|
||||||
|
# XL compiler : consider using -qarch=auto
|
||||||
|
CXXOPT="-O3 -qarch=auto $CXXOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"ifc" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CXXOPT="-O3 $CXXOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"pg" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CXXCOPT="-fast $CXXOPT"
|
||||||
|
# NOTE : PG & Sun use -fast instead -O3
|
||||||
|
elif test "X$psblas_cv_fc" == X"sun" ; then
|
||||||
|
# other compilers ..
|
||||||
|
CXXOPT="-fast $CXXOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"cray" ; then
|
||||||
|
CXXOPT="-O3 $CXXOPT"
|
||||||
|
MPICXX="CC"
|
||||||
|
else
|
||||||
|
CXXOPT="-g -O3 $CXXOPT"
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
# Honor FCFLAGS if they were specified explicitly, but --with-fcopt take precedence
|
||||||
|
if test "X$FCOPT" == "X" ; then
|
||||||
|
FCOPT="$FCFLAGS";
|
||||||
|
fi
|
||||||
|
if test "X$FCOPT" == "X" ; then
|
||||||
|
if test "X$psblas_cv_fc" == "Xgcc" ; then
|
||||||
|
# note that no space should be placed around the equality symbol in assignations
|
||||||
|
# Note : 'native' is valid _only_ on GCC/x86 (32/64 bits)
|
||||||
|
FCOPT="-O3 $FCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"xlf" ; then
|
||||||
|
# XL compiler : consider using -qarch=auto
|
||||||
|
FCOPT="-O3 -qarch=auto -qlanglvl=extended -qxlf2003=polymorphic:autorealloc $FCOPT"
|
||||||
|
FCFLAGS="-qhalt=e -qlanglvl=extended -qxlf2003=polymorphic:autorealloc $FCFLAGS"
|
||||||
|
elif test "X$psblas_cv_fc" == X"ifc" ; then
|
||||||
|
# other compilers ..
|
||||||
|
FCOPT="-O3 $FCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"pg" ; then
|
||||||
|
# other compilers ..
|
||||||
|
FCOPT="-fast $FCOPT"
|
||||||
|
# NOTE : PG & Sun use -fast instead -O3
|
||||||
|
elif test "X$psblas_cv_fc" == X"sun" ; then
|
||||||
|
# other compilers ..
|
||||||
|
FCOPT="-fast $FCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"cray" ; then
|
||||||
|
FCOPT="-O3 -em $FCOPT"
|
||||||
|
elif test "X$psblas_cv_fc" == X"nag" ; then
|
||||||
|
# NAG compiler ..
|
||||||
|
FCOPT="-O2 "
|
||||||
|
# NOTE : PG & Sun use -fast instead -O3
|
||||||
|
else
|
||||||
|
FCOPT="-O2 $FCOPT"
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
if test "X$psblas_cv_fc" == X"nag" ; then
|
||||||
|
# Add needed options
|
||||||
|
FCOPT="$FCOPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
|
||||||
|
EXTRA_OPT="-mismatch_all"
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
# COPT,FCOPT are aliases for CFLAGS,FCFLAGS .
|
||||||
|
|
||||||
|
##############################################################################
|
||||||
|
# Compilers variables selection
|
||||||
|
##############################################################################
|
||||||
|
FC=${FC}
|
||||||
|
CC=${CC}
|
||||||
|
CXX=${CXX}
|
||||||
|
CCOPT="$CCOPT $C99OPT"
|
||||||
|
|
||||||
|
|
||||||
|
##############################################################################
|
||||||
|
# Choice of our compilers, needed by Make.inc
|
||||||
|
##############################################################################
|
||||||
|
if test "X$psblas_cv_fc" == X"cray"
|
||||||
|
then
|
||||||
|
MODEXT=".mod"
|
||||||
|
FMFLAG="-I"
|
||||||
|
FIFLAG="-I"
|
||||||
|
BASEMODNAME=PSB_BASE_MOD
|
||||||
|
PRECMODNAME=PSB_PREC_MOD
|
||||||
|
METHDMODNAME=PSB_KRYLOV_MOD
|
||||||
|
UTILMODNAME=PSB_UTIL_MOD
|
||||||
|
|
||||||
|
else
|
||||||
|
AX_F90_MODULE_EXTENSION
|
||||||
|
AX_F90_MODULE_FLAG
|
||||||
|
MODEXT=".$ax_cv_f90_modext"
|
||||||
|
FMFLAG="${ax_cv_f90_modflag%%[ ]*}"
|
||||||
|
FIFLAG=-I
|
||||||
|
BASEMODNAME=psb_base_mod
|
||||||
|
PRECMODNAME=psb_prec_mod
|
||||||
|
METHDMODNAME=psb_krylov_mod
|
||||||
|
UTILMODNAME=psb_util_mod
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
##############################################################################
|
||||||
|
# Choice of our compilers, needed by Make.inc
|
||||||
|
##############################################################################
|
||||||
|
if test "X$FLINK" == "X" ; then
|
||||||
|
FLINK=${MPF90}
|
||||||
|
fi
|
||||||
|
# Custom test : do we have a module or include for MPI Fortran interface?
|
||||||
|
if test x"$pac_cv_serial_mpi" == x"yes" ; then
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DSERIAL_MPI $psblas_cv_define_prepend-DMPI_MOD $FDEFINES";
|
||||||
|
else
|
||||||
|
PAC_FORTRAN_CHECK_HAVE_MPI_MOD_F08()
|
||||||
|
if test x"$pac_cv_mpi_f08" == x"yes" ; then
|
||||||
|
dnl FDEFINES="$psblas_cv_define_prepend-DMPI_MOD_F08 $FDEFINES";
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DMPI_MOD $FDEFINES";
|
||||||
|
else
|
||||||
|
PAC_FORTRAN_CHECK_HAVE_MPI_MOD(
|
||||||
|
[FDEFINES="$psblas_cv_define_prepend-DMPI_MOD $FDEFINES"],
|
||||||
|
[FDEFINES="$psblas_cv_define_prepend-DMPI_H $FDEFINES"])
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
|
||||||
|
FLINK="$MPIFC"
|
||||||
|
PAC_ARG_OPENMP()
|
||||||
|
if test x"$pac_cv_openmp" == x"yes" ; then
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DOPENMP $FDEFINES";
|
||||||
|
CDEFINES="-DOPENMP $CDEFINES";
|
||||||
|
FCOPT="$FCOPT $pac_cv_openmp_fcopt";
|
||||||
|
CCOPT="$CCOPT $pac_cv_openmp_ccopt";
|
||||||
|
FLINK="$FLINK $pac_cv_openmp_fcopt";
|
||||||
|
fi
|
||||||
|
|
||||||
|
PAC_FORTRAN_HAVE_PSBLAS([AC_MSG_RESULT([yes.])],
|
||||||
|
[AC_MSG_ERROR([no. Could not find working version of PSBLAS.])])
|
||||||
|
|
||||||
|
PAC_FORTRAN_PSBLAS_VERSION()
|
||||||
|
|
||||||
|
if test "x$pac_cv_psblas_major" == "xunknown"; then
|
||||||
|
AC_MSG_ERROR([PSBLAS version major "$pac_cv_psblas_major".])
|
||||||
|
fi
|
||||||
|
if test "x$pac_cv_psblas_minor" == "xunknown"; then
|
||||||
|
AC_MSG_ERROR([PSBLAS version minor "$pac_cv_psblas_minor".])
|
||||||
|
fi
|
||||||
|
if test "x$pac_cv_psblas_patchlevel" == "xunknown"; then
|
||||||
|
AC_MSG_ERROR([PSBLAS patchlevel "$pac_cv_psblas_patchlevel".])
|
||||||
|
fi
|
||||||
|
if (( $pac_cv_psblas_major < 3 )) ||
|
||||||
|
( (( $pac_cv_psblas_major == 3 )) && (( $pac_cv_psblas_minor < 7 ))) ; then
|
||||||
|
AC_MSG_ERROR([I need at least PSBLAS version 3.7.])
|
||||||
|
else
|
||||||
|
AC_MSG_NOTICE([Am configuring with PSBLAS version $pac_cv_psblas_major.$pac_cv_psblas_minor.$pac_cv_psblas_patchlevel.])
|
||||||
|
fi
|
||||||
|
|
||||||
|
PAC_ARG_WITH_IPK
|
||||||
|
PAC_ARG_WITH_LPK
|
||||||
|
# Defaults for IPK/LPK
|
||||||
|
if test x"$pac_cv_ipk_size" == x"" ; then
|
||||||
|
pac_cv_ipk_size=4
|
||||||
|
fi
|
||||||
|
if test x"$pac_cv_lpk_size" == x"" ; then
|
||||||
|
pac_cv_lpk_size=8
|
||||||
|
fi
|
||||||
|
# Enforce sensible combination
|
||||||
|
if (( $pac_cv_lpk_size < $pac_cv_ipk_size )); then
|
||||||
|
AC_MSG_NOTICE([[Invalid combination of size specs IPK ${pac_cv_ipk_size} LPK ${pac_cv_lpk_size}. ]]);
|
||||||
|
AC_MSG_NOTICE([[Forcing equal values]])
|
||||||
|
pac_cv_lpk_size=$pac_cv_ipk_size;
|
||||||
|
fi
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DIPK${pac_cv_ipk_size} $FDEFINES";
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DLPK${pac_cv_lpk_size} $FDEFINES";
|
||||||
|
CDEFINES="-DIPK${pac_cv_ipk_size} -DLPK${pac_cv_lpk_size} $CDEFINES";
|
||||||
|
if test x"$pac_cv_lpk_size" == x8"" ; then
|
||||||
|
CXXDEFINES="-DBIT64 $CXXDEFINES";
|
||||||
|
fi
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Parachute rules for ar and ranlib ... (could cause problems)
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
if test "X$AR" == "X" ; then
|
||||||
|
AR="ar"
|
||||||
|
fi
|
||||||
|
|
||||||
|
if test "X$RANLIB" == "X" ; then
|
||||||
|
RANLIB="ranlib"
|
||||||
|
fi
|
||||||
|
|
||||||
|
# This should be portable
|
||||||
|
AR="${AR} -cur"
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# NOTE :
|
||||||
|
# Missing stuff :
|
||||||
|
# In the case the detected fortran compiler is ifort, icc or gcc
|
||||||
|
# should be valid options.
|
||||||
|
# The same for pg (Portland Group compilers).
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
#
|
||||||
|
# Tests for support of various Fortran features; some of them are critical,
|
||||||
|
# some optional
|
||||||
|
#
|
||||||
|
|
||||||
|
#
|
||||||
|
# Critical features
|
||||||
|
#
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_EXTENDS(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for EXTENDS.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_CLASS_TBP(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for CLASS and type bound procedures.
|
||||||
|
Please get a Fortran compiler that supports them, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_SOURCE(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for SOURCE= allocation.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_HAVE_MOVE_ALLOC(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for MOVE_ALLOC.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_ISO_C_BIND(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for ISO_C_BINDING.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_SAME_TYPE(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for SAME_TYPE_AS.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_EXTENDS_TYPE(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for EXTENDS_TYPE_OF.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_MOLD(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for MOLD= allocation.
|
||||||
|
Please get a Fortran compiler that supports it, e.g. GNU Fortran 4.8.])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_VOLATILE(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for VOLATILE])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_ISO_FORTRAN_ENV(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for ISO_FORTRAN_ENV])]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_FINAL(
|
||||||
|
[],
|
||||||
|
[AC_MSG_ERROR([Sorry, cannot build PSBLAS without support for FINAL])]
|
||||||
|
)
|
||||||
|
|
||||||
|
|
||||||
|
#
|
||||||
|
# Optional features
|
||||||
|
#
|
||||||
|
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_GENERICS(
|
||||||
|
[],
|
||||||
|
[FDEFINES="$psblas_cv_define_prepend-DHAVE_BUGGY_GENERICS $FDEFINES"]
|
||||||
|
)
|
||||||
|
|
||||||
|
PAC_FORTRAN_TEST_FLUSH(
|
||||||
|
[FDEFINES="$psblas_cv_define_prepend-DHAVE_FLUSH_STMT $FDEFINES"],
|
||||||
|
)
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# Additional pathname stuff (yes, it is redundant and confusing...)
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
# -I
|
||||||
|
if test x"$INCLUDEPATH" != "x" ; then
|
||||||
|
FINCLUDES="$FINCLUDES $INCLUDEPATH"
|
||||||
|
CINCLUDES="$CINCLUDES $INCLUDEPATH"
|
||||||
|
fi
|
||||||
|
|
||||||
|
# -L
|
||||||
|
if test x"$LIBRARYPATH" != "x" ; then
|
||||||
|
FINCLUDES="$FINCLUDES $LIBRARYPATH"
|
||||||
|
fi
|
||||||
|
|
||||||
|
# -I
|
||||||
|
if test x"$MODULE_PATH" != "x" ; then
|
||||||
|
FINCLUDES="$FINCLUDES $MODULE_PATH"
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# BLAS library presence checks
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
# Note : The libmkl.a (Intel Math Kernel Library) library could be used, too.
|
||||||
|
# It is sufficient to specify it as -lmkl in the CLIBS or FLIBS or LIBS
|
||||||
|
# and specify its path adjusting -L/path in CFLAGS.
|
||||||
|
|
||||||
|
# Right now it is a matter of user's taste when linking custom applications.
|
||||||
|
# But PSBLAS examples could take advantage of these libraries, too.
|
||||||
|
AC_LANG([Fortran])
|
||||||
|
###############################################################################
|
||||||
|
# BLAS library presence checks
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
# Note : The libmkl.a (Intel Math Kernel Library) library could be used, too.
|
||||||
|
# It is sufficient to specify it as -lmkl in the CLIBS or FLIBS or LIBS
|
||||||
|
# and specify its path adjusting -L/path in CFLAGS.
|
||||||
|
|
||||||
|
# Right now it is a matter of user's taste when linking custom applications.
|
||||||
|
# But PSBLAS examples could take advantage of these libraries, too.
|
||||||
|
|
||||||
|
PAC_BLAS([], [AC_MSG_ERROR([[Cannot find BLAS library, specify a path using --with-blas=DIR/LIB (for example --with-blas=/usr/path/lib/libcxml.a)]])])
|
||||||
|
PAC_LAPACK(
|
||||||
|
[FDEFINES="$psblas_cv_define_prepend-DHAVE_LAPACK $FDEFINES"],
|
||||||
|
)
|
||||||
|
AC_LANG([C])
|
||||||
|
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# BLACS library presence checks
|
||||||
|
###############################################################################
|
||||||
|
#AC_LANG([C])
|
||||||
|
#if test x"$pac_cv_serial_mpi" == x"no" ; then
|
||||||
|
#save_FC="$FC";
|
||||||
|
#save_CC="$CC";
|
||||||
|
#FC="$MPIFC";
|
||||||
|
#CC="$MPICC";
|
||||||
|
#PAC_CHECK_BLACS
|
||||||
|
#FC="$save_FC";
|
||||||
|
#CC="$save_CC";
|
||||||
|
#fi
|
||||||
|
|
||||||
|
PAC_MAKE_IS_GNUMAKE
|
||||||
|
###############################################################################
|
||||||
|
# Auxiliary packages
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
PAC_CHECK_METIS
|
||||||
|
|
||||||
|
AC_MSG_CHECKING([Compatibility between metis and LPK])
|
||||||
|
if test "x$pac_cv_lpk_size" == "x4" ; then
|
||||||
|
if test "x$pac_cv_metis_idx" == "x64" ; then
|
||||||
|
dnl mismatch between metis size and PSBLAS LPK
|
||||||
|
psblas_cv_have_metis="no";
|
||||||
|
dnl
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
if test "x$pac_cv_lpk_size" == "x8" ; then
|
||||||
|
if test "x$pac_cv_metis_idx" == "x32" ; then
|
||||||
|
dnl mismatch between metis size and PSBLAS LPK
|
||||||
|
psblas_cv_have_metis="no";
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
AC_MSG_RESULT([$psblas_cv_have_metis])
|
||||||
|
|
||||||
|
if test "x$pac_cv_metis_idx" == "xunknown" ; then
|
||||||
|
dnl mismatch between metis size and PSBLAS LPK
|
||||||
|
AC_MSG_NOTICE([Unknown METIS bitsize.])
|
||||||
|
$psblas_cv_have_metis = "no";
|
||||||
|
fi
|
||||||
|
if test "x$pac_cv_metis_real" == "xunknown" ; then
|
||||||
|
dnl mismatch between metis size and PSBLAS LPK
|
||||||
|
AC_MSG_NOTICE([Unknown METIS REAL bitsize.])
|
||||||
|
$psblas_cv_have_metis = "no";
|
||||||
|
fi
|
||||||
|
if test "x$psblas_cv_have_metis" == "xyes" ; then
|
||||||
|
FDEFINES="$psblas_cv_define_prepend-DHAVE_METIS $psblas_cv_define_prepend-DMETIS_$pac_cv_metis_idx $psblas_cv_define_prepend-DMETIS_REAL_$pac_cv_metis_real $FDEFINES"
|
||||||
|
CDEFINES="-DHAVE_METIS_ $psblas_cv_metis_includes $CDEFINES -DMETIS_$pac_cv_metis_idx -DMETIS_REAL_$pac_cv_metis_real"
|
||||||
|
METISINCFILE=$psblas_cv_metisincfile
|
||||||
|
fi
|
||||||
|
|
||||||
|
|
||||||
|
PAC_CHECK_MUMPS
|
||||||
|
#
|
||||||
|
# 1. Enable even with LPK=8, internally it will check if
|
||||||
|
# the problem size fits into 4 bytes, very likely since we
|
||||||
|
# are mostly using MUMPS at coarse level.
|
||||||
|
#
|
||||||
|
dnl if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
|
||||||
|
dnl if test "x$pac_cv_psblas_ipk" == "x8" ; then
|
||||||
|
dnl AC_MSG_NOTICE([PSBLAS defines PSB_IPK_ as $pac_cv_psblas_ipk. MUMPS interfacing disabled. ])
|
||||||
|
dnl MUMPS_FLAGS="";
|
||||||
|
dnl MUMPS_LIBS="";
|
||||||
|
dnl amg4psblas_cv_have_mumps=no;
|
||||||
|
dnl fi
|
||||||
|
dnl fi
|
||||||
|
if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
|
||||||
|
if test "x$pac_cv_psblas_lpk" == "x8" ; then
|
||||||
|
AC_MSG_NOTICE([PSBLAS defines PSB_LPK_ as $pac_cv_psblas_lpk. MUMPS interfacing will fail when called in global mode on very large matrices. ])
|
||||||
|
fi
|
||||||
|
if test "x$pac_mumps_fmods_ok" == "xyes" ; then
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DHAVE_MUMPS_ $amg_cv_define_prepend-DHAVE_MUMPS_MODULES_ $MUMPS_MODULES $FDEFINES"
|
||||||
|
MUMPS_FLAGS="-DHave_MUMPS_ $MUMPS_MODULES"
|
||||||
|
elif test "x$pac_mumps_fincs_ok" == "xyes" ; then
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DHAVE_MUMPS_ $amg_cv_define_prepend-DHAVE_MUMPS_INCLUDES_ $MUMPS_FINCLUDES $FDEFINES"
|
||||||
|
MUMPS_FLAGS="-DHave_MUMPS_ $MUMPS_INCLUDES"
|
||||||
|
else
|
||||||
|
# This should not happen
|
||||||
|
MUMPS_FLAGS=""
|
||||||
|
MUMPS_LIBS=""
|
||||||
|
fi
|
||||||
|
else
|
||||||
|
MUMPS_FLAGS=""
|
||||||
|
MUMPS_LIBS=""
|
||||||
|
fi
|
||||||
|
|
||||||
|
PAC_CHECK_UMFPACK
|
||||||
|
if test "x$amg4psblas_cv_have_umfpack" == "xyes" ; then
|
||||||
|
UMF_FLAGS="-DHave_UMF_ $UMF_INCLUDES"
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DHAVE_UMF_ $FDEFINES"
|
||||||
|
else
|
||||||
|
UMF_FLAGS=""
|
||||||
|
fi
|
||||||
|
|
||||||
|
PAC_CHECK_SUPERLU
|
||||||
|
if test "x$amg4psblas_cv_have_superlu" == "xyes" ; then
|
||||||
|
SLU_FLAGS="-DHave_SLU_ -DSLU_VERSION_$pac_slu_version $SLU_INCLUDES"
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DHAVE_SLU_ $FDEFINES"
|
||||||
|
else
|
||||||
|
SLU_FLAGS=""
|
||||||
|
fi
|
||||||
|
|
||||||
|
PAC_CHECK_SUPERLUDIST()
|
||||||
|
|
||||||
|
if test "x$amg4psblas_cv_have_superludist" == "xyes" ; then
|
||||||
|
pac_sludist_version="$amg4psblas_cv_superludist_major";
|
||||||
|
if (($amg4psblas_cv_superludist_major==6)); then
|
||||||
|
if (($amg4psblas_cv_superludist_minor>=3)); then
|
||||||
|
pac_sludist_version="63";
|
||||||
|
fi
|
||||||
|
fi
|
||||||
|
SLUDIST_FLAGS=""
|
||||||
|
SLUDIST_FLAGS="-DHave_SLUDist_ -DSLUD_VERSION_$pac_sludist_version $SLUDIST_INCLUDES"
|
||||||
|
FDEFINES="$amg_cv_define_prepend-DHAVE_SLUDIST_ $FDEFINES"
|
||||||
|
else
|
||||||
|
SLUDIST_FLAGS=""
|
||||||
|
fi
|
||||||
|
|
||||||
|
##############################################
|
||||||
|
FINCLUDES="$PSBLAS_INCLUDES"
|
||||||
|
|
||||||
|
AMGFDEFINES="$FDEFINES"
|
||||||
|
AMGCDEFINES="$CDEFINES"
|
||||||
|
AMGCXXDEFINES="$CXXDEFINES"
|
||||||
|
|
||||||
|
LIBDIR=lib
|
||||||
|
BASELIBNAME=libpsb_base.a
|
||||||
|
PRECLIBNAME=libpsb_prec.a
|
||||||
|
METHDLIBNAME=libpsb_krylov.a
|
||||||
|
UTILLIBNAME=libpsb_util.a
|
||||||
|
AMGLIBNAME=libamg_prec.a
|
||||||
|
|
||||||
|
COMPILERULES='
|
||||||
|
PSBLDLIBS=$(LAPACK) $(BLAS) $(METIS_LIB) $(AMD_LIB) $(LIBS)
|
||||||
|
CXXDEFINES=$(PSBCXXDEFINES)
|
||||||
|
CDEFINES=$(PSBCDEFINES)
|
||||||
|
FDEFINES=$(PSBFDEFINES)
|
||||||
|
|
||||||
|
|
||||||
|
# These should be portable rules, arent they?
|
||||||
|
.c.o:
|
||||||
|
$(CC) $(CCOPT) $(CINCLUDES) $(CDEFINES) -c $< -o $@
|
||||||
|
.f90.o:
|
||||||
|
$(FC) $(FCOPT) $(FINCLUDES) -c $< -o $@
|
||||||
|
.F90.o:
|
||||||
|
$(FC) $(FCOPT) $(FINCLUDES) $(FDEFINES) -c $< -o $@
|
||||||
|
.cpp.o:
|
||||||
|
$(CXX) $(CXXOPT) $(CXXINCLUDES) $(CXXDEFINES) -c $< -o $@'
|
||||||
|
###############################################################################
|
||||||
|
# Variable substitutions : the Make.inc.in will have these @VARIABLES@
|
||||||
|
# substituted.
|
||||||
|
|
||||||
|
AC_SUBST(PSBLAS_DIR)
|
||||||
|
AC_SUBST(PSBLAS_INCDIR)
|
||||||
|
AC_SUBST(PSBLAS_MODDIR)
|
||||||
|
AC_SUBST(PSBLAS_LIBDIR)
|
||||||
|
AC_SUBST(PSBLAS_INCLUDES)
|
||||||
|
dnl AC_SUBST(PSBLAS_INSTALL_MAKEINC)
|
||||||
|
AC_SUBST(PSBLAS_LIBS)
|
||||||
|
AC_SUBST(PSBLAS_RULES)
|
||||||
|
|
||||||
|
AC_SUBST(INSTALL)
|
||||||
|
AC_SUBST(INSTALL_DATA)
|
||||||
|
AC_SUBST(INSTALL_DIR)
|
||||||
|
AC_SUBST(INSTALL_LIBDIR)
|
||||||
|
AC_SUBST(INSTALL_INCLUDEDIR)
|
||||||
|
AC_SUBST(INSTALL_MODULESDIR)
|
||||||
|
AC_SUBST(INSTALL_DOCSDIR)
|
||||||
|
AC_SUBST(INSTALL_SAMPLESDIR)
|
||||||
|
|
||||||
|
AC_SUBST(EXTRA_LIBS)
|
||||||
|
AC_SUBST(BLAS_LIBS)
|
||||||
|
AC_SUBST(LAPACK_LIBS)
|
||||||
|
AC_SUBST(METIS_LIBS)
|
||||||
|
AC_SUBST(MUMPS_FLAGS)
|
||||||
|
AC_SUBST(MUMPS_LIBS)
|
||||||
|
AC_SUBST(SLU_FLAGS)
|
||||||
|
AC_SUBST(SLU_LIBS)
|
||||||
|
AC_SUBST(UMF_FLAGS)
|
||||||
|
AC_SUBST(UMF_LIBS)
|
||||||
|
AC_SUBST(SLUDIST_FLAGS)
|
||||||
|
AC_SUBST(SLUDIST_LIBS)
|
||||||
|
AC_SUBST(AMGFDEFINES)
|
||||||
|
AC_SUBST(AMGCDEFINES)
|
||||||
|
AC_SUBST(AMGCXXDEFINES)
|
||||||
|
AC_SUBST(MODEXT)
|
||||||
|
AC_SUBST(COMPILERULES)
|
||||||
|
AC_SUBST(FDEFINES)
|
||||||
|
AC_SUBST(CDEFINES)
|
||||||
|
AC_SUBST(BASEMODNAME)
|
||||||
|
AC_SUBST(PRECMODNAME)
|
||||||
|
AC_SUBST(METHDMODNAME)
|
||||||
|
AC_SUBST(UTILMODNAME)
|
||||||
|
AC_SUBST(AMGLIBNAME)
|
||||||
|
AC_SUBST(MPIFC)
|
||||||
|
AC_SUBST(MPICC)
|
||||||
|
AC_SUBST(MPICXX)
|
||||||
|
AC_SUBST(FCOPT)
|
||||||
|
AC_SUBST(CCOPT)
|
||||||
|
AC_SUBST(CXXOPT)
|
||||||
|
AC_SUBST(EXTRA_OPT)
|
||||||
|
AC_SUBST(FAKEMPI)
|
||||||
|
AC_SUBST(FIFLAG)
|
||||||
|
AC_SUBST(FMFLAG)
|
||||||
|
AC_SUBST(MODEXT)
|
||||||
|
AC_SUBST(FLINK)
|
||||||
|
AC_SUBST(LIBS)
|
||||||
|
AC_SUBST(AR)
|
||||||
|
AC_SUBST(RANLIB)
|
||||||
|
AC_SUBST(MPIFC)
|
||||||
|
AC_SUBST(MPIFCC)
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
# the following files will be created by Automake
|
||||||
|
|
||||||
|
AC_CONFIG_FILES([Make_n.inc])
|
||||||
|
AC_OUTPUT()
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
dnl Please note that brackets around variable identifiers are absolutely needed for compatibility..
|
||||||
|
AC_MSG_NOTICE([
|
||||||
|
${PACKAGE_NAME} ${amg4psblas_cv_version} has been configured as follows:
|
||||||
|
|
||||||
|
PSBLAS library : ${PSBLAS_DIR}
|
||||||
|
MUMPS detected : ${amg4psblas_cv_have_mumps}
|
||||||
|
SuperLU detected : ${amg4psblas_cv_have_superlu}
|
||||||
|
SuperLU_Dist detected : ${amg4psblas_cv_have_superludist}
|
||||||
|
UMFPack detected : ${amg4psblas_cv_have_umfpack}
|
||||||
|
|
||||||
|
If you are satisfied, run 'make' to build ${PACKAGE_NAME} and its documentation; otherwise
|
||||||
|
type ./configure --help=short for a complete list of configure options specific to ${PACKAGE_NAME}.
|
||||||
|
dnl To install the program and its documentation, run 'make install' if you are root,
|
||||||
|
dnl or run 'su -c "make install"' if you are not root.
|
||||||
|
])
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
@ -0,0 +1,931 @@
|
|||||||
|
module amg_d_genpde_mod
|
||||||
|
|
||||||
|
|
||||||
|
use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,&
|
||||||
|
& psb_dspmat_type, psb_d_vect_type, dzero, done,&
|
||||||
|
& psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type
|
||||||
|
|
||||||
|
interface
|
||||||
|
function d_func_3d(x,y,z) result(val)
|
||||||
|
import :: psb_dpk_
|
||||||
|
real(psb_dpk_), intent(in) :: x,y,z
|
||||||
|
real(psb_dpk_) :: val
|
||||||
|
end function d_func_3d
|
||||||
|
end interface
|
||||||
|
|
||||||
|
interface amg_gen_pde3d
|
||||||
|
module procedure amg_d_gen_pde3d
|
||||||
|
end interface amg_gen_pde3d
|
||||||
|
|
||||||
|
interface
|
||||||
|
function d_func_2d(x,y) result(val)
|
||||||
|
import :: psb_dpk_
|
||||||
|
real(psb_dpk_), intent(in) :: x,y
|
||||||
|
real(psb_dpk_) :: val
|
||||||
|
end function d_func_2d
|
||||||
|
end interface
|
||||||
|
|
||||||
|
interface amg_gen_pde2d
|
||||||
|
module procedure amg_d_gen_pde2d
|
||||||
|
end interface amg_gen_pde2d
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
function d_null_func_2d(x,y) result(val)
|
||||||
|
|
||||||
|
real(psb_dpk_), intent(in) :: x,y
|
||||||
|
real(psb_dpk_) :: val
|
||||||
|
|
||||||
|
val = dzero
|
||||||
|
|
||||||
|
end function d_null_func_2d
|
||||||
|
|
||||||
|
function d_null_func_3d(x,y,z) result(val)
|
||||||
|
|
||||||
|
real(psb_dpk_), intent(in) :: x,y,z
|
||||||
|
real(psb_dpk_) :: val
|
||||||
|
|
||||||
|
val = dzero
|
||||||
|
|
||||||
|
end function d_null_func_3d
|
||||||
|
|
||||||
|
!
|
||||||
|
! subroutine to allocate and fill in the coefficient matrix and
|
||||||
|
! the rhs.
|
||||||
|
!
|
||||||
|
subroutine amg_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,&
|
||||||
|
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv)
|
||||||
|
use psb_base_mod
|
||||||
|
use psb_util_mod
|
||||||
|
!
|
||||||
|
! Discretizes the partial differential equation
|
||||||
|
!
|
||||||
|
! d a1 d(u) d a1 d(u) d a1 d(u) b1 d(u) b2 d(u) b3 d(u)
|
||||||
|
! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f
|
||||||
|
! dx dx dy dy dz dz dx dy dz
|
||||||
|
!
|
||||||
|
! with Dirichlet boundary conditions
|
||||||
|
! u = g
|
||||||
|
!
|
||||||
|
! on the unit cube 0<=x,y,z<=1.
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g
|
||||||
|
integer(psb_ipk_) :: idim
|
||||||
|
type(psb_dspmat_type) :: a
|
||||||
|
type(psb_d_vect_type) :: xv,bv
|
||||||
|
type(psb_desc_type) :: desc_a
|
||||||
|
integer(psb_ipk_) :: info
|
||||||
|
type(psb_ctxt_type) :: ctxt
|
||||||
|
character :: afmt*5
|
||||||
|
procedure(d_func_3d), optional :: f
|
||||||
|
class(psb_d_base_sparse_mat), optional :: amold
|
||||||
|
class(psb_d_base_vect_type), optional :: vmold
|
||||||
|
integer(psb_ipk_), optional :: partition, nrl,iv(:)
|
||||||
|
|
||||||
|
! Local variables.
|
||||||
|
|
||||||
|
integer(psb_ipk_), parameter :: nb=20
|
||||||
|
type(psb_d_csc_sparse_mat) :: acsc
|
||||||
|
type(psb_d_coo_sparse_mat) :: acoo
|
||||||
|
type(psb_d_csr_sparse_mat) :: acsr
|
||||||
|
real(psb_dpk_) :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
|
||||||
|
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
|
||||||
|
integer(psb_lpk_) :: m,n,glob_row,nt
|
||||||
|
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
|
||||||
|
! For 3D partition
|
||||||
|
! Note: integer control variables going directly into an MPI call
|
||||||
|
! must be 4 bytes, i.e. psb_mpk_
|
||||||
|
integer(psb_mpk_) :: npdims(3), npp, minfo
|
||||||
|
integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz
|
||||||
|
integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
|
||||||
|
! Process grid
|
||||||
|
integer(psb_ipk_) :: np, iam
|
||||||
|
integer(psb_ipk_) :: icoeff
|
||||||
|
integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:)
|
||||||
|
real(psb_dpk_), allocatable :: val(:)
|
||||||
|
! deltah dimension of each grid cell
|
||||||
|
! deltat discretization time
|
||||||
|
real(psb_dpk_) :: deltah, sqdeltah, deltah2
|
||||||
|
real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero
|
||||||
|
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
|
||||||
|
integer(psb_ipk_) :: err_act
|
||||||
|
procedure(d_func_3d), pointer :: f_
|
||||||
|
character(len=20) :: name, ch_err,tmpfmt
|
||||||
|
|
||||||
|
info = psb_success_
|
||||||
|
name = 'd_create_matrix'
|
||||||
|
call psb_erractionsave(err_act)
|
||||||
|
|
||||||
|
call psb_info(ctxt, iam, np)
|
||||||
|
|
||||||
|
|
||||||
|
if (present(f)) then
|
||||||
|
f_ => f
|
||||||
|
else
|
||||||
|
f_ => d_null_func_3d
|
||||||
|
end if
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
deltah = done/(idim+2)
|
||||||
|
sqdeltah = deltah*deltah
|
||||||
|
deltah2 = 2.0_psb_dpk_* deltah
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
|
||||||
|
! initialize array descriptor and sparse matrix storage. provide an
|
||||||
|
! estimate of the number of non zeroes
|
||||||
|
|
||||||
|
m = (1_psb_lpk_*idim)*idim*idim
|
||||||
|
n = m
|
||||||
|
nnz = 7*((n+np-1)/np)
|
||||||
|
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
|
||||||
|
t0 = psb_wtime()
|
||||||
|
select case(partition_)
|
||||||
|
case(1)
|
||||||
|
! A BLOCK partition
|
||||||
|
if (present(nrl)) then
|
||||||
|
nr = nrl
|
||||||
|
else
|
||||||
|
!
|
||||||
|
! Using a simple BLOCK distribution.
|
||||||
|
!
|
||||||
|
nt = (m+np-1)/np
|
||||||
|
nr = max(0,min(nt,m-(iam*nt)))
|
||||||
|
end if
|
||||||
|
|
||||||
|
nt = nr
|
||||||
|
call psb_sum(ctxt,nt)
|
||||||
|
if (nt /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! First example of use of CDALL: specify for each process a number of
|
||||||
|
! contiguous rows
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,nl=nr)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
! A partition defined by the user through IV
|
||||||
|
|
||||||
|
if (present(iv)) then
|
||||||
|
if (size(iv) /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Second example of use of CDALL: specify for each row the
|
||||||
|
! process that owns it
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vg=iv)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(3)
|
||||||
|
! A 3-dimensional partition
|
||||||
|
|
||||||
|
! A nifty MPI function will split the process list
|
||||||
|
npdims = 0
|
||||||
|
call mpi_dims_create(np,3,npdims,info)
|
||||||
|
npx = npdims(1)
|
||||||
|
npy = npdims(2)
|
||||||
|
npz = npdims(3)
|
||||||
|
|
||||||
|
allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
|
||||||
|
! We can reuse idx2ijk for process indices as well.
|
||||||
|
call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
|
||||||
|
! Now let's split the 3D cube in hexahedra
|
||||||
|
call dist1Didx(bndx,idim,npx)
|
||||||
|
mynx = bndx(iamx+1)-bndx(iamx)
|
||||||
|
call dist1Didx(bndy,idim,npy)
|
||||||
|
myny = bndy(iamy+1)-bndy(iamy)
|
||||||
|
call dist1Didx(bndz,idim,npz)
|
||||||
|
mynz = bndz(iamz+1)-bndz(iamz)
|
||||||
|
|
||||||
|
! How many indices do I own?
|
||||||
|
nlr = mynx*myny*mynz
|
||||||
|
allocate(myidx(nlr))
|
||||||
|
! Now, let's generate the list of indices I own
|
||||||
|
nr = 0
|
||||||
|
do i=bndx(iamx),bndx(iamx+1)-1
|
||||||
|
do j=bndy(iamy),bndy(iamy+1)-1
|
||||||
|
do k=bndz(iamz),bndz(iamz+1)-1
|
||||||
|
nr = nr + 1
|
||||||
|
call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
if (nr /= nlr) then
|
||||||
|
write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
|
||||||
|
& nr,nlr,mynx,myny,mynz
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Third example of use of CDALL: specify for each process
|
||||||
|
! the set of global indices it owns.
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vl=myidx)
|
||||||
|
|
||||||
|
!
|
||||||
|
! Specify process topology
|
||||||
|
!
|
||||||
|
block
|
||||||
|
!
|
||||||
|
! Use adjcncy methods
|
||||||
|
!
|
||||||
|
integer(psb_mpk_), allocatable :: neighbours(:)
|
||||||
|
integer(psb_mpk_) :: cnt
|
||||||
|
logical, parameter :: debug_adj=.true.
|
||||||
|
if (debug_adj.and.(np > 1)) then
|
||||||
|
cnt = 0
|
||||||
|
allocate(neighbours(np))
|
||||||
|
if (iamx < npx-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx+1,iamy,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy < npy-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy+1,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamz < npz-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy,iamz+1,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamx >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx-1,iamy,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy-1,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamz >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy,iamz-1,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
call psb_realloc(cnt, neighbours,info)
|
||||||
|
call desc_a%set_p_adjcncy(neighbours)
|
||||||
|
!write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
|
||||||
|
end if
|
||||||
|
end block
|
||||||
|
|
||||||
|
case default
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end select
|
||||||
|
|
||||||
|
|
||||||
|
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
|
||||||
|
! define rhs from boundary conditions; also build initial guess
|
||||||
|
if (info == psb_success_) call psb_geall(xv,desc_a,info)
|
||||||
|
if (info == psb_success_) call psb_geall(bv,desc_a,info)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
talc = psb_wtime()-t0
|
||||||
|
|
||||||
|
if (info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='allocation rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
! we build an auxiliary matrix consisting of one row at a
|
||||||
|
! time; just a small matrix. might be extended to generate
|
||||||
|
! a bunch of rows per call.
|
||||||
|
!
|
||||||
|
allocate(val(20*nb),irow(20*nb),&
|
||||||
|
&icol(20*nb),stat=info)
|
||||||
|
if (info /= psb_success_ ) then
|
||||||
|
info=psb_err_alloc_dealloc_
|
||||||
|
call psb_errpush(info,name)
|
||||||
|
goto 9999
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
! loop over rows belonging to current process in a block
|
||||||
|
! distribution.
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
do ii=1, nlr,nb
|
||||||
|
ib = min(nb,nlr-ii+1)
|
||||||
|
icoeff = 1
|
||||||
|
do k=1,ib
|
||||||
|
i=ii+k-1
|
||||||
|
! local matrix pointer
|
||||||
|
glob_row=myidx(i)
|
||||||
|
! compute gridpoint coordinates
|
||||||
|
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
|
||||||
|
! x, y, z coordinates
|
||||||
|
x = (ix-1)*deltah
|
||||||
|
y = (iy-1)*deltah
|
||||||
|
z = (iz-1)*deltah
|
||||||
|
zt(k) = f_(x,y,z)
|
||||||
|
! internal point: build discretization
|
||||||
|
!
|
||||||
|
! term depending on (x-1,y,z)
|
||||||
|
!
|
||||||
|
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
|
||||||
|
if (ix == 1) then
|
||||||
|
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y-1,z)
|
||||||
|
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
|
||||||
|
if (iy == 1) then
|
||||||
|
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y,z-1)
|
||||||
|
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
|
||||||
|
if (iz == 1) then
|
||||||
|
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
! term depending on (x,y,z)
|
||||||
|
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
|
||||||
|
& + c(x,y,z)
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
! term depending on (x,y,z+1)
|
||||||
|
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
|
||||||
|
if (iz == idim) then
|
||||||
|
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y+1,z)
|
||||||
|
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
|
||||||
|
if (iy == idim) then
|
||||||
|
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x+1,y,z)
|
||||||
|
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
|
||||||
|
if (ix==idim) then
|
||||||
|
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
end do
|
||||||
|
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
zt(:)=dzero
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
end do
|
||||||
|
|
||||||
|
tgen = psb_wtime()-t1
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='insert rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
deallocate(val,irow,icol)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
call psb_cdasb(desc_a,info)
|
||||||
|
tcdasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
if (info == psb_success_) then
|
||||||
|
if (present(amold)) then
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold)
|
||||||
|
else
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
|
||||||
|
if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
tasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
ttot = psb_wtime() - t0
|
||||||
|
|
||||||
|
call psb_amx(ctxt,talc)
|
||||||
|
call psb_amx(ctxt,tgen)
|
||||||
|
call psb_amx(ctxt,tasb)
|
||||||
|
call psb_amx(ctxt,ttot)
|
||||||
|
if(iam == psb_root_) then
|
||||||
|
tmpfmt = a%get_fmt()
|
||||||
|
write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||||
|
& tmpfmt
|
||||||
|
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
|
||||||
|
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
|
||||||
|
write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb
|
||||||
|
write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb
|
||||||
|
write(psb_out_unit,'("-total time : ",es12.5)') ttot
|
||||||
|
|
||||||
|
end if
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
return
|
||||||
|
|
||||||
|
9999 call psb_error_handler(ctxt,err_act)
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine amg_d_gen_pde3d
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
!
|
||||||
|
! subroutine to allocate and fill in the coefficient matrix and
|
||||||
|
! the rhs.
|
||||||
|
!
|
||||||
|
subroutine amg_d_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,&
|
||||||
|
& a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv)
|
||||||
|
use psb_base_mod
|
||||||
|
use psb_util_mod
|
||||||
|
!
|
||||||
|
! Discretizes the partial differential equation
|
||||||
|
!
|
||||||
|
! d d(u) d d(u) b1 d(u) b2 d(u)
|
||||||
|
! - -- a1 ---- - -- a1 ---- + ----- + ------ + c u = f
|
||||||
|
! dx dx dy dy dx dy
|
||||||
|
!
|
||||||
|
! with Dirichlet boundary conditions
|
||||||
|
! u = g
|
||||||
|
!
|
||||||
|
! on the unit square 0<=x,y<=1.
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
procedure(d_func_2d) :: b1,b2,c,a1,a2,g
|
||||||
|
integer(psb_ipk_) :: idim
|
||||||
|
type(psb_dspmat_type) :: a
|
||||||
|
type(psb_d_vect_type) :: xv,bv
|
||||||
|
type(psb_desc_type) :: desc_a
|
||||||
|
integer(psb_ipk_) :: info
|
||||||
|
type(psb_ctxt_type) :: ctxt
|
||||||
|
character :: afmt*5
|
||||||
|
procedure(d_func_2d), optional :: f
|
||||||
|
class(psb_d_base_sparse_mat), optional :: amold
|
||||||
|
class(psb_d_base_vect_type), optional :: vmold
|
||||||
|
integer(psb_ipk_), optional :: partition, nrl,iv(:)
|
||||||
|
! Local variables.
|
||||||
|
|
||||||
|
integer(psb_ipk_), parameter :: nb=20
|
||||||
|
type(psb_d_csc_sparse_mat) :: acsc
|
||||||
|
type(psb_d_coo_sparse_mat) :: acoo
|
||||||
|
type(psb_d_csr_sparse_mat) :: acsr
|
||||||
|
real(psb_dpk_) :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
|
||||||
|
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
|
||||||
|
integer(psb_lpk_) :: m,n,glob_row,nt
|
||||||
|
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
|
||||||
|
! For 2D partition
|
||||||
|
! Note: integer control variables going directly into an MPI call
|
||||||
|
! must be 4 bytes, i.e. psb_mpk_
|
||||||
|
integer(psb_mpk_) :: npdims(2), npp, minfo
|
||||||
|
integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny
|
||||||
|
integer(psb_ipk_), allocatable :: bndx(:),bndy(:)
|
||||||
|
! Process grid
|
||||||
|
integer(psb_ipk_) :: np, iam
|
||||||
|
integer(psb_ipk_) :: icoeff
|
||||||
|
integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:)
|
||||||
|
real(psb_dpk_), allocatable :: val(:)
|
||||||
|
! deltah dimension of each grid cell
|
||||||
|
! deltat discretization time
|
||||||
|
real(psb_dpk_) :: deltah, sqdeltah, deltah2, dd
|
||||||
|
real(psb_dpk_), parameter :: rhs=0.d0,one=done,zero=0.d0
|
||||||
|
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
|
||||||
|
integer(psb_ipk_) :: err_act
|
||||||
|
procedure(d_func_2d), pointer :: f_
|
||||||
|
character(len=20) :: name, ch_err,tmpfmt
|
||||||
|
|
||||||
|
info = psb_success_
|
||||||
|
name = 'create_matrix'
|
||||||
|
call psb_erractionsave(err_act)
|
||||||
|
|
||||||
|
call psb_info(ctxt, iam, np)
|
||||||
|
|
||||||
|
|
||||||
|
if (present(f)) then
|
||||||
|
f_ => f
|
||||||
|
else
|
||||||
|
f_ => d_null_func_2d
|
||||||
|
end if
|
||||||
|
|
||||||
|
deltah = done/(idim+2)
|
||||||
|
sqdeltah = deltah*deltah
|
||||||
|
deltah2 = 2.0_psb_dpk_* deltah
|
||||||
|
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
|
||||||
|
! initialize array descriptor and sparse matrix storage. provide an
|
||||||
|
! estimate of the number of non zeroes
|
||||||
|
|
||||||
|
m = (1_psb_lpk_)*idim*idim
|
||||||
|
n = m
|
||||||
|
nnz = 7*((n+np-1)/np)
|
||||||
|
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
|
||||||
|
t0 = psb_wtime()
|
||||||
|
select case(partition_)
|
||||||
|
case(1)
|
||||||
|
! A BLOCK partition
|
||||||
|
if (present(nrl)) then
|
||||||
|
nr = nrl
|
||||||
|
else
|
||||||
|
!
|
||||||
|
! Using a simple BLOCK distribution.
|
||||||
|
!
|
||||||
|
nt = (m+np-1)/np
|
||||||
|
nr = max(0,min(nt,m-(iam*nt)))
|
||||||
|
end if
|
||||||
|
|
||||||
|
nt = nr
|
||||||
|
call psb_sum(ctxt,nt)
|
||||||
|
if (nt /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! First example of use of CDALL: specify for each process a number of
|
||||||
|
! contiguous rows
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,nl=nr)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
! A partition defined by the user through IV
|
||||||
|
|
||||||
|
if (present(iv)) then
|
||||||
|
if (size(iv) /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Second example of use of CDALL: specify for each row the
|
||||||
|
! process that owns it
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vg=iv)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(3)
|
||||||
|
! A 2-dimensional partition
|
||||||
|
|
||||||
|
! A nifty MPI function will split the process list
|
||||||
|
npdims = 0
|
||||||
|
call mpi_dims_create(np,2,npdims,info)
|
||||||
|
npx = npdims(1)
|
||||||
|
npy = npdims(2)
|
||||||
|
|
||||||
|
allocate(bndx(0:npx),bndy(0:npy))
|
||||||
|
! We can reuse idx2ijk for process indices as well.
|
||||||
|
call idx2ijk(iamx,iamy,iam,npx,npy,base=0)
|
||||||
|
! Now let's split the 2D square in rectangles
|
||||||
|
call dist1Didx(bndx,idim,npx)
|
||||||
|
mynx = bndx(iamx+1)-bndx(iamx)
|
||||||
|
call dist1Didx(bndy,idim,npy)
|
||||||
|
myny = bndy(iamy+1)-bndy(iamy)
|
||||||
|
|
||||||
|
! How many indices do I own?
|
||||||
|
nlr = mynx*myny
|
||||||
|
allocate(myidx(nlr))
|
||||||
|
! Now, let's generate the list of indices I own
|
||||||
|
nr = 0
|
||||||
|
do i=bndx(iamx),bndx(iamx+1)-1
|
||||||
|
do j=bndy(iamy),bndy(iamy+1)-1
|
||||||
|
nr = nr + 1
|
||||||
|
call ijk2idx(myidx(nr),i,j,idim,idim)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
if (nr /= nlr) then
|
||||||
|
write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',&
|
||||||
|
& nr,nlr,mynx,myny
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Third example of use of CDALL: specify for each process
|
||||||
|
! the set of global indices it owns.
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vl=myidx)
|
||||||
|
|
||||||
|
!
|
||||||
|
! Specify process topology
|
||||||
|
!
|
||||||
|
block
|
||||||
|
!
|
||||||
|
! Use adjcncy methods
|
||||||
|
!
|
||||||
|
integer(psb_mpk_), allocatable :: neighbours(:)
|
||||||
|
integer(psb_mpk_) :: cnt
|
||||||
|
logical, parameter :: debug_adj=.true.
|
||||||
|
if (debug_adj.and.(np > 1)) then
|
||||||
|
cnt = 0
|
||||||
|
allocate(neighbours(np))
|
||||||
|
if (iamx < npx-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx+1,iamy,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy < npy-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy+1,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamx >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx-1,iamy,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy-1,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
call psb_realloc(cnt, neighbours,info)
|
||||||
|
call desc_a%set_p_adjcncy(neighbours)
|
||||||
|
!write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
|
||||||
|
end if
|
||||||
|
end block
|
||||||
|
|
||||||
|
case default
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end select
|
||||||
|
|
||||||
|
|
||||||
|
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
|
||||||
|
! define rhs from boundary conditions; also build initial guess
|
||||||
|
if (info == psb_success_) call psb_geall(xv,desc_a,info)
|
||||||
|
if (info == psb_success_) call psb_geall(bv,desc_a,info)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
talc = psb_wtime()-t0
|
||||||
|
|
||||||
|
if (info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='allocation rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
! we build an auxiliary matrix consisting of one row at a
|
||||||
|
! time; just a small matrix. might be extended to generate
|
||||||
|
! a bunch of rows per call.
|
||||||
|
!
|
||||||
|
allocate(val(20*nb),irow(20*nb),&
|
||||||
|
&icol(20*nb),stat=info)
|
||||||
|
if (info /= psb_success_ ) then
|
||||||
|
info=psb_err_alloc_dealloc_
|
||||||
|
call psb_errpush(info,name)
|
||||||
|
goto 9999
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
! loop over rows belonging to current process in a block
|
||||||
|
! distribution.
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
do ii=1, nlr,nb
|
||||||
|
ib = min(nb,nlr-ii+1)
|
||||||
|
icoeff = 1
|
||||||
|
do k=1,ib
|
||||||
|
i=ii+k-1
|
||||||
|
! local matrix pointer
|
||||||
|
glob_row=myidx(i)
|
||||||
|
! compute gridpoint coordinates
|
||||||
|
call idx2ijk(ix,iy,glob_row,idim,idim)
|
||||||
|
! x, y coordinates
|
||||||
|
x = (ix-1)*deltah
|
||||||
|
y = (iy-1)*deltah
|
||||||
|
|
||||||
|
zt(k) = f_(x,y)
|
||||||
|
! internal point: build discretization
|
||||||
|
!
|
||||||
|
! term depending on (x-1,y)
|
||||||
|
!
|
||||||
|
val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
|
||||||
|
if (ix == 1) then
|
||||||
|
zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y-1)
|
||||||
|
val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
|
||||||
|
if (iy == 1) then
|
||||||
|
zt(k) = g(x,dzero)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
! term depending on (x,y)
|
||||||
|
val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
! term depending on (x,y+1)
|
||||||
|
val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
|
||||||
|
if (iy == idim) then
|
||||||
|
zt(k) = g(x,done)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x+1,y)
|
||||||
|
val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
|
||||||
|
if (ix==idim) then
|
||||||
|
zt(k) = g(done,y)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
end do
|
||||||
|
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
zt(:)=dzero
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
end do
|
||||||
|
|
||||||
|
tgen = psb_wtime()-t1
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='insert rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
deallocate(val,irow,icol)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
call psb_cdasb(desc_a,info)
|
||||||
|
tcdasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
if (info == psb_success_) then
|
||||||
|
if (present(amold)) then
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold)
|
||||||
|
else
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
|
||||||
|
if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
tasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
ttot = psb_wtime() - t0
|
||||||
|
|
||||||
|
call psb_amx(ctxt,talc)
|
||||||
|
call psb_amx(ctxt,tgen)
|
||||||
|
call psb_amx(ctxt,tasb)
|
||||||
|
call psb_amx(ctxt,ttot)
|
||||||
|
if(iam == psb_root_) then
|
||||||
|
tmpfmt = a%get_fmt()
|
||||||
|
write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||||
|
& tmpfmt
|
||||||
|
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
|
||||||
|
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
|
||||||
|
write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb
|
||||||
|
write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb
|
||||||
|
write(psb_out_unit,'("-total time : ",es12.5)') ttot
|
||||||
|
|
||||||
|
end if
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
return
|
||||||
|
|
||||||
|
9999 continue
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
if (err_act == psb_act_abort_) then
|
||||||
|
call psb_error(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
return
|
||||||
|
end subroutine amg_d_gen_pde2d
|
||||||
|
end module amg_d_genpde_mod
|
||||||
@ -0,0 +1,931 @@
|
|||||||
|
module amg_s_genpde_mod
|
||||||
|
|
||||||
|
|
||||||
|
use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,&
|
||||||
|
& psb_sspmat_type, psb_s_vect_type, szero, sone,&
|
||||||
|
& psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type
|
||||||
|
|
||||||
|
interface
|
||||||
|
function s_func_3d(x,y,z) result(val)
|
||||||
|
import :: psb_spk_
|
||||||
|
real(psb_spk_), intent(in) :: x,y,z
|
||||||
|
real(psb_spk_) :: val
|
||||||
|
end function s_func_3d
|
||||||
|
end interface
|
||||||
|
|
||||||
|
interface amg_gen_pde3d
|
||||||
|
module procedure amg_s_gen_pde3d
|
||||||
|
end interface amg_gen_pde3d
|
||||||
|
|
||||||
|
interface
|
||||||
|
function s_func_2d(x,y) result(val)
|
||||||
|
import :: psb_spk_
|
||||||
|
real(psb_spk_), intent(in) :: x,y
|
||||||
|
real(psb_spk_) :: val
|
||||||
|
end function s_func_2d
|
||||||
|
end interface
|
||||||
|
|
||||||
|
interface amg_gen_pde2d
|
||||||
|
module procedure amg_s_gen_pde2d
|
||||||
|
end interface amg_gen_pde2d
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
function s_null_func_2d(x,y) result(val)
|
||||||
|
|
||||||
|
real(psb_spk_), intent(in) :: x,y
|
||||||
|
real(psb_spk_) :: val
|
||||||
|
|
||||||
|
val = szero
|
||||||
|
|
||||||
|
end function s_null_func_2d
|
||||||
|
|
||||||
|
function s_null_func_3d(x,y,z) result(val)
|
||||||
|
|
||||||
|
real(psb_spk_), intent(in) :: x,y,z
|
||||||
|
real(psb_spk_) :: val
|
||||||
|
|
||||||
|
val = szero
|
||||||
|
|
||||||
|
end function s_null_func_3d
|
||||||
|
|
||||||
|
!
|
||||||
|
! subroutine to allocate and fill in the coefficient matrix and
|
||||||
|
! the rhs.
|
||||||
|
!
|
||||||
|
subroutine amg_s_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,&
|
||||||
|
& a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv)
|
||||||
|
use psb_base_mod
|
||||||
|
use psb_util_mod
|
||||||
|
!
|
||||||
|
! Discretizes the partial differential equation
|
||||||
|
!
|
||||||
|
! d a1 d(u) d a1 d(u) d a1 d(u) b1 d(u) b2 d(u) b3 d(u)
|
||||||
|
! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f
|
||||||
|
! dx dx dy dy dz dz dx dy dz
|
||||||
|
!
|
||||||
|
! with Dirichlet boundary conditions
|
||||||
|
! u = g
|
||||||
|
!
|
||||||
|
! on the unit cube 0<=x,y,z<=1.
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
procedure(s_func_3d) :: b1,b2,b3,c,a1,a2,a3,g
|
||||||
|
integer(psb_ipk_) :: idim
|
||||||
|
type(psb_sspmat_type) :: a
|
||||||
|
type(psb_s_vect_type) :: xv,bv
|
||||||
|
type(psb_desc_type) :: desc_a
|
||||||
|
integer(psb_ipk_) :: info
|
||||||
|
type(psb_ctxt_type) :: ctxt
|
||||||
|
character :: afmt*5
|
||||||
|
procedure(s_func_3d), optional :: f
|
||||||
|
class(psb_s_base_sparse_mat), optional :: amold
|
||||||
|
class(psb_s_base_vect_type), optional :: vmold
|
||||||
|
integer(psb_ipk_), optional :: partition, nrl,iv(:)
|
||||||
|
|
||||||
|
! Local variables.
|
||||||
|
|
||||||
|
integer(psb_ipk_), parameter :: nb=20
|
||||||
|
type(psb_s_csc_sparse_mat) :: acsc
|
||||||
|
type(psb_s_coo_sparse_mat) :: acoo
|
||||||
|
type(psb_s_csr_sparse_mat) :: acsr
|
||||||
|
real(psb_spk_) :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
|
||||||
|
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
|
||||||
|
integer(psb_lpk_) :: m,n,glob_row,nt
|
||||||
|
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
|
||||||
|
! For 3D partition
|
||||||
|
! Note: integer control variables going directly into an MPI call
|
||||||
|
! must be 4 bytes, i.e. psb_mpk_
|
||||||
|
integer(psb_mpk_) :: npdims(3), npp, minfo
|
||||||
|
integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz
|
||||||
|
integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
|
||||||
|
! Process grid
|
||||||
|
integer(psb_ipk_) :: np, iam
|
||||||
|
integer(psb_ipk_) :: icoeff
|
||||||
|
integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:)
|
||||||
|
real(psb_spk_), allocatable :: val(:)
|
||||||
|
! deltah dimension of each grid cell
|
||||||
|
! deltat discretization time
|
||||||
|
real(psb_spk_) :: deltah, sqdeltah, deltah2
|
||||||
|
real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero
|
||||||
|
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
|
||||||
|
integer(psb_ipk_) :: err_act
|
||||||
|
procedure(s_func_3d), pointer :: f_
|
||||||
|
character(len=20) :: name, ch_err,tmpfmt
|
||||||
|
|
||||||
|
info = psb_success_
|
||||||
|
name = 's_create_matrix'
|
||||||
|
call psb_erractionsave(err_act)
|
||||||
|
|
||||||
|
call psb_info(ctxt, iam, np)
|
||||||
|
|
||||||
|
|
||||||
|
if (present(f)) then
|
||||||
|
f_ => f
|
||||||
|
else
|
||||||
|
f_ => s_null_func_3d
|
||||||
|
end if
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
deltah = sone/(idim+2)
|
||||||
|
sqdeltah = deltah*deltah
|
||||||
|
deltah2 = 2.0_psb_spk_* deltah
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
|
||||||
|
! initialize array descriptor and sparse matrix storage. provide an
|
||||||
|
! estimate of the number of non zeroes
|
||||||
|
|
||||||
|
m = (1_psb_lpk_*idim)*idim*idim
|
||||||
|
n = m
|
||||||
|
nnz = 7*((n+np-1)/np)
|
||||||
|
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
|
||||||
|
t0 = psb_wtime()
|
||||||
|
select case(partition_)
|
||||||
|
case(1)
|
||||||
|
! A BLOCK partition
|
||||||
|
if (present(nrl)) then
|
||||||
|
nr = nrl
|
||||||
|
else
|
||||||
|
!
|
||||||
|
! Using a simple BLOCK distribution.
|
||||||
|
!
|
||||||
|
nt = (m+np-1)/np
|
||||||
|
nr = max(0,min(nt,m-(iam*nt)))
|
||||||
|
end if
|
||||||
|
|
||||||
|
nt = nr
|
||||||
|
call psb_sum(ctxt,nt)
|
||||||
|
if (nt /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! First example of use of CDALL: specify for each process a number of
|
||||||
|
! contiguous rows
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,nl=nr)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
! A partition defined by the user through IV
|
||||||
|
|
||||||
|
if (present(iv)) then
|
||||||
|
if (size(iv) /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Second example of use of CDALL: specify for each row the
|
||||||
|
! process that owns it
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vg=iv)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(3)
|
||||||
|
! A 3-dimensional partition
|
||||||
|
|
||||||
|
! A nifty MPI function will split the process list
|
||||||
|
npdims = 0
|
||||||
|
call mpi_dims_create(np,3,npdims,info)
|
||||||
|
npx = npdims(1)
|
||||||
|
npy = npdims(2)
|
||||||
|
npz = npdims(3)
|
||||||
|
|
||||||
|
allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
|
||||||
|
! We can reuse idx2ijk for process indices as well.
|
||||||
|
call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
|
||||||
|
! Now let's split the 3D cube in hexahedra
|
||||||
|
call dist1Didx(bndx,idim,npx)
|
||||||
|
mynx = bndx(iamx+1)-bndx(iamx)
|
||||||
|
call dist1Didx(bndy,idim,npy)
|
||||||
|
myny = bndy(iamy+1)-bndy(iamy)
|
||||||
|
call dist1Didx(bndz,idim,npz)
|
||||||
|
mynz = bndz(iamz+1)-bndz(iamz)
|
||||||
|
|
||||||
|
! How many indices do I own?
|
||||||
|
nlr = mynx*myny*mynz
|
||||||
|
allocate(myidx(nlr))
|
||||||
|
! Now, let's generate the list of indices I own
|
||||||
|
nr = 0
|
||||||
|
do i=bndx(iamx),bndx(iamx+1)-1
|
||||||
|
do j=bndy(iamy),bndy(iamy+1)-1
|
||||||
|
do k=bndz(iamz),bndz(iamz+1)-1
|
||||||
|
nr = nr + 1
|
||||||
|
call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
if (nr /= nlr) then
|
||||||
|
write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
|
||||||
|
& nr,nlr,mynx,myny,mynz
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Third example of use of CDALL: specify for each process
|
||||||
|
! the set of global indices it owns.
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vl=myidx)
|
||||||
|
|
||||||
|
!
|
||||||
|
! Specify process topology
|
||||||
|
!
|
||||||
|
block
|
||||||
|
!
|
||||||
|
! Use adjcncy methods
|
||||||
|
!
|
||||||
|
integer(psb_mpk_), allocatable :: neighbours(:)
|
||||||
|
integer(psb_mpk_) :: cnt
|
||||||
|
logical, parameter :: debug_adj=.true.
|
||||||
|
if (debug_adj.and.(np > 1)) then
|
||||||
|
cnt = 0
|
||||||
|
allocate(neighbours(np))
|
||||||
|
if (iamx < npx-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx+1,iamy,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy < npy-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy+1,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamz < npz-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy,iamz+1,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamx >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx-1,iamy,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy-1,iamz,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
if (iamz >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy,iamz-1,npx,npy,npz,base=0)
|
||||||
|
end if
|
||||||
|
call psb_realloc(cnt, neighbours,info)
|
||||||
|
call desc_a%set_p_adjcncy(neighbours)
|
||||||
|
!write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
|
||||||
|
end if
|
||||||
|
end block
|
||||||
|
|
||||||
|
case default
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end select
|
||||||
|
|
||||||
|
|
||||||
|
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
|
||||||
|
! define rhs from boundary conditions; also build initial guess
|
||||||
|
if (info == psb_success_) call psb_geall(xv,desc_a,info)
|
||||||
|
if (info == psb_success_) call psb_geall(bv,desc_a,info)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
talc = psb_wtime()-t0
|
||||||
|
|
||||||
|
if (info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='allocation rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
! we build an auxiliary matrix consisting of one row at a
|
||||||
|
! time; just a small matrix. might be extended to generate
|
||||||
|
! a bunch of rows per call.
|
||||||
|
!
|
||||||
|
allocate(val(20*nb),irow(20*nb),&
|
||||||
|
&icol(20*nb),stat=info)
|
||||||
|
if (info /= psb_success_ ) then
|
||||||
|
info=psb_err_alloc_dealloc_
|
||||||
|
call psb_errpush(info,name)
|
||||||
|
goto 9999
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
! loop over rows belonging to current process in a block
|
||||||
|
! distribution.
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
do ii=1, nlr,nb
|
||||||
|
ib = min(nb,nlr-ii+1)
|
||||||
|
icoeff = 1
|
||||||
|
do k=1,ib
|
||||||
|
i=ii+k-1
|
||||||
|
! local matrix pointer
|
||||||
|
glob_row=myidx(i)
|
||||||
|
! compute gridpoint coordinates
|
||||||
|
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
|
||||||
|
! x, y, z coordinates
|
||||||
|
x = (ix-1)*deltah
|
||||||
|
y = (iy-1)*deltah
|
||||||
|
z = (iz-1)*deltah
|
||||||
|
zt(k) = f_(x,y,z)
|
||||||
|
! internal point: build discretization
|
||||||
|
!
|
||||||
|
! term depending on (x-1,y,z)
|
||||||
|
!
|
||||||
|
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
|
||||||
|
if (ix == 1) then
|
||||||
|
zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y-1,z)
|
||||||
|
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
|
||||||
|
if (iy == 1) then
|
||||||
|
zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y,z-1)
|
||||||
|
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
|
||||||
|
if (iz == 1) then
|
||||||
|
zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
! term depending on (x,y,z)
|
||||||
|
val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
|
||||||
|
& + c(x,y,z)
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
! term depending on (x,y,z+1)
|
||||||
|
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
|
||||||
|
if (iz == idim) then
|
||||||
|
zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y+1,z)
|
||||||
|
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
|
||||||
|
if (iy == idim) then
|
||||||
|
zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x+1,y,z)
|
||||||
|
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
|
||||||
|
if (ix==idim) then
|
||||||
|
zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
end do
|
||||||
|
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
zt(:)=szero
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
end do
|
||||||
|
|
||||||
|
tgen = psb_wtime()-t1
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='insert rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
deallocate(val,irow,icol)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
call psb_cdasb(desc_a,info)
|
||||||
|
tcdasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
if (info == psb_success_) then
|
||||||
|
if (present(amold)) then
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold)
|
||||||
|
else
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
|
||||||
|
if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
tasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
ttot = psb_wtime() - t0
|
||||||
|
|
||||||
|
call psb_amx(ctxt,talc)
|
||||||
|
call psb_amx(ctxt,tgen)
|
||||||
|
call psb_amx(ctxt,tasb)
|
||||||
|
call psb_amx(ctxt,ttot)
|
||||||
|
if(iam == psb_root_) then
|
||||||
|
tmpfmt = a%get_fmt()
|
||||||
|
write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||||
|
& tmpfmt
|
||||||
|
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
|
||||||
|
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
|
||||||
|
write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb
|
||||||
|
write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb
|
||||||
|
write(psb_out_unit,'("-total time : ",es12.5)') ttot
|
||||||
|
|
||||||
|
end if
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
return
|
||||||
|
|
||||||
|
9999 call psb_error_handler(ctxt,err_act)
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine amg_s_gen_pde3d
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
!
|
||||||
|
! subroutine to allocate and fill in the coefficient matrix and
|
||||||
|
! the rhs.
|
||||||
|
!
|
||||||
|
subroutine amg_s_gen_pde2d(ctxt,idim,a,bv,xv,desc_a,afmt,&
|
||||||
|
& a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv)
|
||||||
|
use psb_base_mod
|
||||||
|
use psb_util_mod
|
||||||
|
!
|
||||||
|
! Discretizes the partial differential equation
|
||||||
|
!
|
||||||
|
! d d(u) d d(u) b1 d(u) b2 d(u)
|
||||||
|
! - -- a1 ---- - -- a1 ---- + ----- + ------ + c u = f
|
||||||
|
! dx dx dy dy dx dy
|
||||||
|
!
|
||||||
|
! with Dirichlet boundary conditions
|
||||||
|
! u = g
|
||||||
|
!
|
||||||
|
! on the unit square 0<=x,y<=1.
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! Note that if b1=b2=c=0., the PDE is the Laplace equation.
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
procedure(s_func_2d) :: b1,b2,c,a1,a2,g
|
||||||
|
integer(psb_ipk_) :: idim
|
||||||
|
type(psb_sspmat_type) :: a
|
||||||
|
type(psb_s_vect_type) :: xv,bv
|
||||||
|
type(psb_desc_type) :: desc_a
|
||||||
|
integer(psb_ipk_) :: info
|
||||||
|
type(psb_ctxt_type) :: ctxt
|
||||||
|
character :: afmt*5
|
||||||
|
procedure(s_func_2d), optional :: f
|
||||||
|
class(psb_s_base_sparse_mat), optional :: amold
|
||||||
|
class(psb_s_base_vect_type), optional :: vmold
|
||||||
|
integer(psb_ipk_), optional :: partition, nrl,iv(:)
|
||||||
|
! Local variables.
|
||||||
|
|
||||||
|
integer(psb_ipk_), parameter :: nb=20
|
||||||
|
type(psb_s_csc_sparse_mat) :: acsc
|
||||||
|
type(psb_s_coo_sparse_mat) :: acoo
|
||||||
|
type(psb_s_csr_sparse_mat) :: acsr
|
||||||
|
real(psb_spk_) :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
|
||||||
|
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
|
||||||
|
integer(psb_lpk_) :: m,n,glob_row,nt
|
||||||
|
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
|
||||||
|
! For 2D partition
|
||||||
|
! Note: integer control variables going directly into an MPI call
|
||||||
|
! must be 4 bytes, i.e. psb_mpk_
|
||||||
|
integer(psb_mpk_) :: npdims(2), npp, minfo
|
||||||
|
integer(psb_ipk_) :: npx,npy,iamx,iamy,mynx,myny
|
||||||
|
integer(psb_ipk_), allocatable :: bndx(:),bndy(:)
|
||||||
|
! Process grid
|
||||||
|
integer(psb_ipk_) :: np, iam
|
||||||
|
integer(psb_ipk_) :: icoeff
|
||||||
|
integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:)
|
||||||
|
real(psb_spk_), allocatable :: val(:)
|
||||||
|
! deltah dimension of each grid cell
|
||||||
|
! deltat discretization time
|
||||||
|
real(psb_spk_) :: deltah, sqdeltah, deltah2, dd
|
||||||
|
real(psb_spk_), parameter :: rhs=0.d0,one=sone,zero=0.d0
|
||||||
|
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
|
||||||
|
integer(psb_ipk_) :: err_act
|
||||||
|
procedure(s_func_2d), pointer :: f_
|
||||||
|
character(len=20) :: name, ch_err,tmpfmt
|
||||||
|
|
||||||
|
info = psb_success_
|
||||||
|
name = 'create_matrix'
|
||||||
|
call psb_erractionsave(err_act)
|
||||||
|
|
||||||
|
call psb_info(ctxt, iam, np)
|
||||||
|
|
||||||
|
|
||||||
|
if (present(f)) then
|
||||||
|
f_ => f
|
||||||
|
else
|
||||||
|
f_ => s_null_func_2d
|
||||||
|
end if
|
||||||
|
|
||||||
|
deltah = sone/(idim+2)
|
||||||
|
sqdeltah = deltah*deltah
|
||||||
|
deltah2 = 2.0_psb_spk_* deltah
|
||||||
|
|
||||||
|
|
||||||
|
if (present(partition)) then
|
||||||
|
if ((1<= partition).and.(partition <= 3)) then
|
||||||
|
partition_ = partition
|
||||||
|
else
|
||||||
|
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
partition_ = 3
|
||||||
|
end if
|
||||||
|
|
||||||
|
! initialize array descriptor and sparse matrix storage. provide an
|
||||||
|
! estimate of the number of non zeroes
|
||||||
|
|
||||||
|
m = (1_psb_lpk_)*idim*idim
|
||||||
|
n = m
|
||||||
|
nnz = 7*((n+np-1)/np)
|
||||||
|
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
|
||||||
|
t0 = psb_wtime()
|
||||||
|
select case(partition_)
|
||||||
|
case(1)
|
||||||
|
! A BLOCK partition
|
||||||
|
if (present(nrl)) then
|
||||||
|
nr = nrl
|
||||||
|
else
|
||||||
|
!
|
||||||
|
! Using a simple BLOCK distribution.
|
||||||
|
!
|
||||||
|
nt = (m+np-1)/np
|
||||||
|
nr = max(0,min(nt,m-(iam*nt)))
|
||||||
|
end if
|
||||||
|
|
||||||
|
nt = nr
|
||||||
|
call psb_sum(ctxt,nt)
|
||||||
|
if (nt /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! First example of use of CDALL: specify for each process a number of
|
||||||
|
! contiguous rows
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,nl=nr)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(2)
|
||||||
|
! A partition defined by the user through IV
|
||||||
|
|
||||||
|
if (present(iv)) then
|
||||||
|
if (size(iv) /= m) then
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Second example of use of CDALL: specify for each row the
|
||||||
|
! process that owns it
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vg=iv)
|
||||||
|
myidx = desc_a%get_global_indices()
|
||||||
|
nlr = size(myidx)
|
||||||
|
|
||||||
|
case(3)
|
||||||
|
! A 2-dimensional partition
|
||||||
|
|
||||||
|
! A nifty MPI function will split the process list
|
||||||
|
npdims = 0
|
||||||
|
call mpi_dims_create(np,2,npdims,info)
|
||||||
|
npx = npdims(1)
|
||||||
|
npy = npdims(2)
|
||||||
|
|
||||||
|
allocate(bndx(0:npx),bndy(0:npy))
|
||||||
|
! We can reuse idx2ijk for process indices as well.
|
||||||
|
call idx2ijk(iamx,iamy,iam,npx,npy,base=0)
|
||||||
|
! Now let's split the 2D square in rectangles
|
||||||
|
call dist1Didx(bndx,idim,npx)
|
||||||
|
mynx = bndx(iamx+1)-bndx(iamx)
|
||||||
|
call dist1Didx(bndy,idim,npy)
|
||||||
|
myny = bndy(iamy+1)-bndy(iamy)
|
||||||
|
|
||||||
|
! How many indices do I own?
|
||||||
|
nlr = mynx*myny
|
||||||
|
allocate(myidx(nlr))
|
||||||
|
! Now, let's generate the list of indices I own
|
||||||
|
nr = 0
|
||||||
|
do i=bndx(iamx),bndx(iamx+1)-1
|
||||||
|
do j=bndy(iamy),bndy(iamy+1)-1
|
||||||
|
nr = nr + 1
|
||||||
|
call ijk2idx(myidx(nr),i,j,idim,idim)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
if (nr /= nlr) then
|
||||||
|
write(psb_err_unit,*) iam,iamx,iamy, 'Initialization error: NR vs NLR ',&
|
||||||
|
& nr,nlr,mynx,myny
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
end if
|
||||||
|
|
||||||
|
!
|
||||||
|
! Third example of use of CDALL: specify for each process
|
||||||
|
! the set of global indices it owns.
|
||||||
|
!
|
||||||
|
call psb_cdall(ctxt,desc_a,info,vl=myidx)
|
||||||
|
|
||||||
|
!
|
||||||
|
! Specify process topology
|
||||||
|
!
|
||||||
|
block
|
||||||
|
!
|
||||||
|
! Use adjcncy methods
|
||||||
|
!
|
||||||
|
integer(psb_mpk_), allocatable :: neighbours(:)
|
||||||
|
integer(psb_mpk_) :: cnt
|
||||||
|
logical, parameter :: debug_adj=.true.
|
||||||
|
if (debug_adj.and.(np > 1)) then
|
||||||
|
cnt = 0
|
||||||
|
allocate(neighbours(np))
|
||||||
|
if (iamx < npx-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx+1,iamy,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy < npy-1) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy+1,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamx >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx-1,iamy,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
if (iamy >0) then
|
||||||
|
cnt = cnt + 1
|
||||||
|
call ijk2idx(neighbours(cnt),iamx,iamy-1,npx,npy,base=0)
|
||||||
|
end if
|
||||||
|
call psb_realloc(cnt, neighbours,info)
|
||||||
|
call desc_a%set_p_adjcncy(neighbours)
|
||||||
|
!write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
|
||||||
|
end if
|
||||||
|
end block
|
||||||
|
|
||||||
|
case default
|
||||||
|
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
|
||||||
|
info = -1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
call psb_abort(ctxt)
|
||||||
|
return
|
||||||
|
end select
|
||||||
|
|
||||||
|
|
||||||
|
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
|
||||||
|
! define rhs from boundary conditions; also build initial guess
|
||||||
|
if (info == psb_success_) call psb_geall(xv,desc_a,info)
|
||||||
|
if (info == psb_success_) call psb_geall(bv,desc_a,info)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
talc = psb_wtime()-t0
|
||||||
|
|
||||||
|
if (info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='allocation rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
! we build an auxiliary matrix consisting of one row at a
|
||||||
|
! time; just a small matrix. might be extended to generate
|
||||||
|
! a bunch of rows per call.
|
||||||
|
!
|
||||||
|
allocate(val(20*nb),irow(20*nb),&
|
||||||
|
&icol(20*nb),stat=info)
|
||||||
|
if (info /= psb_success_ ) then
|
||||||
|
info=psb_err_alloc_dealloc_
|
||||||
|
call psb_errpush(info,name)
|
||||||
|
goto 9999
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
! loop over rows belonging to current process in a block
|
||||||
|
! distribution.
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
do ii=1, nlr,nb
|
||||||
|
ib = min(nb,nlr-ii+1)
|
||||||
|
icoeff = 1
|
||||||
|
do k=1,ib
|
||||||
|
i=ii+k-1
|
||||||
|
! local matrix pointer
|
||||||
|
glob_row=myidx(i)
|
||||||
|
! compute gridpoint coordinates
|
||||||
|
call idx2ijk(ix,iy,glob_row,idim,idim)
|
||||||
|
! x, y coordinates
|
||||||
|
x = (ix-1)*deltah
|
||||||
|
y = (iy-1)*deltah
|
||||||
|
|
||||||
|
zt(k) = f_(x,y)
|
||||||
|
! internal point: build discretization
|
||||||
|
!
|
||||||
|
! term depending on (x-1,y)
|
||||||
|
!
|
||||||
|
val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
|
||||||
|
if (ix == 1) then
|
||||||
|
zt(k) = g(szero,y)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x,y-1)
|
||||||
|
val(icoeff) = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
|
||||||
|
if (iy == 1) then
|
||||||
|
zt(k) = g(x,szero)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
! term depending on (x,y)
|
||||||
|
val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
! term depending on (x,y+1)
|
||||||
|
val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
|
||||||
|
if (iy == idim) then
|
||||||
|
zt(k) = g(x,sone)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
! term depending on (x+1,y)
|
||||||
|
val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
|
||||||
|
if (ix==idim) then
|
||||||
|
zt(k) = g(sone,y)*(-val(icoeff)) + zt(k)
|
||||||
|
else
|
||||||
|
call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
|
||||||
|
irow(icoeff) = glob_row
|
||||||
|
icoeff = icoeff+1
|
||||||
|
endif
|
||||||
|
|
||||||
|
end do
|
||||||
|
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
zt(:)=szero
|
||||||
|
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
|
||||||
|
if(info /= psb_success_) exit
|
||||||
|
end do
|
||||||
|
|
||||||
|
tgen = psb_wtime()-t1
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='insert rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
|
||||||
|
deallocate(val,irow,icol)
|
||||||
|
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
call psb_cdasb(desc_a,info)
|
||||||
|
tcdasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
t1 = psb_wtime()
|
||||||
|
if (info == psb_success_) then
|
||||||
|
if (present(amold)) then
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold)
|
||||||
|
else
|
||||||
|
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
|
||||||
|
if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
|
||||||
|
if(info /= psb_success_) then
|
||||||
|
info=psb_err_from_subroutine_
|
||||||
|
ch_err='asb rout.'
|
||||||
|
call psb_errpush(info,name,a_err=ch_err)
|
||||||
|
goto 9999
|
||||||
|
end if
|
||||||
|
tasb = psb_wtime()-t1
|
||||||
|
call psb_barrier(ctxt)
|
||||||
|
ttot = psb_wtime() - t0
|
||||||
|
|
||||||
|
call psb_amx(ctxt,talc)
|
||||||
|
call psb_amx(ctxt,tgen)
|
||||||
|
call psb_amx(ctxt,tasb)
|
||||||
|
call psb_amx(ctxt,ttot)
|
||||||
|
if(iam == psb_root_) then
|
||||||
|
tmpfmt = a%get_fmt()
|
||||||
|
write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||||
|
& tmpfmt
|
||||||
|
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
|
||||||
|
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
|
||||||
|
write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb
|
||||||
|
write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb
|
||||||
|
write(psb_out_unit,'("-total time : ",es12.5)') ttot
|
||||||
|
|
||||||
|
end if
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
return
|
||||||
|
|
||||||
|
9999 continue
|
||||||
|
call psb_erractionrestore(err_act)
|
||||||
|
if (err_act == psb_act_abort_) then
|
||||||
|
call psb_error(ctxt)
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
return
|
||||||
|
end subroutine amg_s_gen_pde2d
|
||||||
|
end module amg_s_genpde_mod
|
||||||
Loading…
Reference in New Issue