From: Wolfgang Bangerth
Date: Tue, 30 Jan 2001 16:38:28 +0000 (+0000)
Subject: Include support for the Harwell Subroutine Library (HSL) into the library, in particu...
X-Git-Tag: v8.0.0~19756
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9c7852abd1f83263d070564274f001baf4ee5cd;p=dealii.git
Include support for the Harwell Subroutine Library (HSL) into the library, in particular for the two sparse direct solvers MA27 and MA47.
git-svn-id: https://svn.dealii.org/trunk@3844 0785d39b-7218-0410-832d-ea1e28bc413d
---
diff --git a/deal.II/Makefile b/deal.II/Makefile
index 7f863cbed4..2a116b96b2 100644
--- a/deal.II/Makefile
+++ b/deal.II/Makefile
@@ -15,8 +15,9 @@ help:
@echo "= base lac: debug and optimized libraries base and lac ="
@echo "= 1dg 2dg 3dg: (1) and deal.II debug version for specified dimension="
@echo "= 1d 2d 3d: (2) and deal.II for specified dimension ="
- @echo "= all: base, lac, 1d, 2d and 3d ="
- @echo "= debug: baseg, lacg, 1dg, 2dg and 3dg ="
+ @echo "= contrib: make additional libraries in contrib, if there are any ="
+ @echo "= all: base, lac, 1d, 2d and 3d, contrib ="
+ @echo "= debug: baseg, lacg, 1dg, 2dg and 3dg, contrib ="
@echo "= online-doc: generate the documentation in HTML format ="
@echo "= printable-doc: generate the documentation in Postscript format. ="
@echo "= You need to have Doc++ and LaTeX installed for that. ="
@@ -71,7 +72,7 @@ lac:
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 3dg
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 3d
-all:
+all: contrib
cd $D/base ; $(MAKE) $(MAKEOPTIONS) libg
cd $D/base ; $(MAKE) $(MAKEOPTIONS) lib
cd $D/lac ; $(MAKE) $(MAKEOPTIONS) libg
@@ -83,13 +84,16 @@ all:
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 3dg
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 3d
-debug:
+debug: contrib
cd $D/base ; $(MAKE) $(MAKEOPTIONS) libg
cd $D/lac ; $(MAKE) $(MAKEOPTIONS) libg
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 2dg
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 1dg
cd $D/deal.II ; $(MAKE) $(MAKEOPTIONS) 3dg
+contrib:
+ cd $D/contrib ; $(MAKE)
+
online-doc:
cd $D/doc ; $(MAKE)
@@ -163,6 +167,6 @@ clean-examples:
.PHONY: base baseg lac lacg 1d 2d 3d 1dg 2dg 3dg all \
- online-doc printable-doc tex-doc \
+ online-doc printable-doc tex-doc contrib \
clean clean-contrib clean-base clean-lac clean-dealII \
clean-doc clean-examples TAGS TODO
diff --git a/deal.II/common/Make.global_options.in b/deal.II/common/Make.global_options.in
index 911c08df1a..1b5c0d16a3 100644
--- a/deal.II/common/Make.global_options.in
+++ b/deal.II/common/Make.global_options.in
@@ -33,21 +33,25 @@ with-multithreading = @withmultithreading@
enable-shared = @enableshared@
DEAL_II_MAJOR = @DEAL_II_MAJOR@
DEAL_II_MINOR = @DEAL_II_MINOR@
+USE_CONTRIB_HSL = @USE_CONTRIB_HSL@
+# if HSL sublibs are used, add F77 libs to linker flags
+ifeq ($(USE_CONTRIB_HSL),yes)
+ LDFLAGS += $(F77LDFLAGS)
+endif
######################################################
# now configuration for compiler flags, pathes, etc
# path for libraries
-lib-path-base = $D/base/lib
-lib-path-lac = $D/lac/lib
-lib-path-deal2 = $D/deal.II/lib
+lib-path-base = $D/base/lib
+lib-path-lac = $D/lac/lib
+lib-path-deal2 = $D/deal.II/lib
+lib-path-contrib-hsl = $D/contrib/hsl/lib
-LIBPATH += $(addprefix -L,$(lib-path-base) $(lib-path-lac) $(lib-path-deal2))
+LIBPATH += $(addprefix -L,$(lib-path-base) $(lib-path-lac) $(lib-path-deal2))
-# DEAL.II libraries are found by the Makefile, so no -L option is needed
-#LIBS += $(LIBPATH)
# add search path for libraries for make dependencies
@@ -55,6 +59,16 @@ vpath %.a $(lib-path-base)
vpath %.a $(lib-path-lac)
vpath %.a $(lib-path-deal2)
+# if HSL sublibs are used, add them to the path
+ifeq ($(USE_CONTRIB_HSL),yes)
+ LIBPATH += $(addprefix -L,$(lib-path-contrib-hsl))
+ vpath %.a $(lib-path-contrib-hsl)
+endif
+
+
+
+
+
lib-suffix = @lib_suffix@
lib-base.o = $(lib-path-base)/libbase$(lib-suffix)
@@ -71,15 +85,27 @@ lib-deal2-1d.g = $(lib-path-deal2)/libdeal_II_1d.g$(lib-suffix)
lib-deal2-2d.g = $(lib-path-deal2)/libdeal_II_2d.g$(lib-suffix)
lib-deal2-3d.g = $(lib-path-deal2)/libdeal_II_3d.g$(lib-suffix)
+lib-contrib-hsl = $(lib-path-contrib-hsl)/libhsl$(lib-suffix)
+
+# if HSL is used, always link it in conjunction with liblac
+ifeq ($(USE_CONTRIB_HSL),yes)
+ lib-lac.g += $(lib-contrib-hsl)
+ lib-lac.o += $(lib-contrib-hsl)
+endif
+
# include paths. add library paths to existing variable
-include-path-base = $D/base/include
-include-path-lac = $D/lac/include
-include-path-deal2 = $D/deal.II/include
+include-path-base = $D/base/include
+include-path-lac = $D/lac/include
+include-path-deal2 = $D/deal.II/include
+include-path-contrib-hsl = $D/contrib/hsl/include
+
+
INCLUDE += $(addprefix -I, $(include-path-base) \
$(include-path-lac) \
- $(include-path-deal2))
+ $(include-path-deal2)\
+ $(include-path-contrib-hsl))
# compiler flags for debug and optimized mode
CXXFLAGS.g = @DEFS@ @CXXFLAGSG@ $(INCLUDE)
diff --git a/deal.II/configure b/deal.II/configure
index 15a2f2a5c6..31cb7e90a9 100755
--- a/deal.II/configure
+++ b/deal.II/configure
@@ -1925,12 +1925,44 @@ EOF
fi
+
+
+echo $ac_n "checking for HSL subroutines""... $ac_c" 1>&6
+echo "configure:1932: checking for HSL subroutines" >&5
+hsl_subroutines=""
+if test -r contrib/hsl/source/ma27.f ; then
+ hsl_subroutines="$hsl_subroutines MA27"
+ cat >> confdefs.h <<\EOF
+#define HAVE_HSL_MA27 1
+EOF
+
+fi
+
+if test -r contrib/hsl/source/ma47.f -a -r contrib/hsl/source/ma47dep.f ; then
+ hsl_subroutines="$hsl_subroutines MA47"
+ cat >> confdefs.h <<\EOF
+#define HAVE_HSL_MA47 1
+EOF
+
+fi
+
+echo "$ac_t""$hsl_subroutines" 1>&6
+
+if test "x$hsl_subroutines" != "x" ; then
+ USE_CONTRIB_HSL=yes
+else
+ USE_CONTRIB_HSL=no
+fi
+
+
+
+
# Check whether --with-blas or --without-blas was given.
if test "${with_blas+set}" = set; then
withval="$with_blas"
if test "$withval" != no ; then
echo $ac_n "checking for BLAS""... $ac_c" 1>&6
-echo "configure:1934: checking for BLAS" >&5
+echo "configure:1966: checking for BLAS" >&5
if test -r "$withval" ; then
echo "$ac_t""found" 1>&6
else
@@ -1946,7 +1978,7 @@ EOF
else
echo $ac_n "checking for main in -lblas""... $ac_c" 1>&6
-echo "configure:1950: checking for main in -lblas" >&5
+echo "configure:1982: checking for main in -lblas" >&5
ac_lib_var=`echo blas'_'main | sed 'y%./+-%__p_%'`
if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then
echo $ac_n "(cached) $ac_c" 1>&6
@@ -1954,14 +1986,14 @@ else
ac_save_LIBS="$LIBS"
LIBS="-lblas $LIBS"
cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then
+if { (eval echo configure:1997: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then
rm -rf conftest*
eval "ac_cv_lib_$ac_lib_var=yes"
else
@@ -1993,7 +2025,7 @@ if test "${with_lapack+set}" = set; then
withval="$with_lapack"
if test "$withval" != no ; then
echo $ac_n "checking for LAPACK""... $ac_c" 1>&6
-echo "configure:1997: checking for LAPACK" >&5
+echo "configure:2029: checking for LAPACK" >&5
if test -r "$withval" ; then
echo "$ac_t""found" 1>&6
else
@@ -2009,7 +2041,7 @@ EOF
else
echo $ac_n "checking for dgels_ in -llapack""... $ac_c" 1>&6
-echo "configure:2013: checking for dgels_ in -llapack" >&5
+echo "configure:2045: checking for dgels_ in -llapack" >&5
ac_lib_var=`echo lapack'_'dgels_ | sed 'y%./+-%__p_%'`
if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then
echo $ac_n "(cached) $ac_c" 1>&6
@@ -2017,7 +2049,7 @@ else
ac_save_LIBS="$LIBS"
LIBS="-llapack $LIBS"
cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then
+if { (eval echo configure:2067: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then
rm -rf conftest*
eval "ac_cv_lib_$ac_lib_var=yes"
else
@@ -2098,7 +2130,7 @@ else
fi
echo $ac_n "checking for kdoc""... $ac_c" 1>&6
-echo "configure:2102: checking for kdoc" >&5
+echo "configure:2134: checking for kdoc" >&5
if test "$kdocdir" != ${DEAL2_DIR}/contrib/kdoc/bin ; then
if test -r $kdocdir/kdoc ; then
echo "$ac_t""found" 1>&6
@@ -2134,7 +2166,7 @@ if test "$docxx" = to-be-determined ; then
# Extract the first word of ""doc++"", so it can be a program name with args.
set dummy "doc++"; ac_word=$2
echo $ac_n "checking for $ac_word""... $ac_c" 1>&6
-echo "configure:2138: checking for $ac_word" >&5
+echo "configure:2170: checking for $ac_word" >&5
if eval "test \"`echo '$''{'ac_cv_path_docxx'+set}'`\" = set"; then
echo $ac_n "(cached) $ac_c" 1>&6
else
@@ -2168,7 +2200,7 @@ fi
else
echo $ac_n "checking for doc++""... $ac_c" 1>&6
-echo "configure:2172: checking for doc++" >&5
+echo "configure:2204: checking for doc++" >&5
if test -x "$docxx" ; then
echo "$ac_t""yes" 1>&6
else
@@ -2183,7 +2215,7 @@ fi
echo $ac_n "checking whether getrusage is properly declared""... $ac_c" 1>&6
-echo "configure:2187: checking whether getrusage is properly declared" >&5
+echo "configure:2219: checking whether getrusage is properly declared" >&5
ac_ext=C
# CXXFLAGS is not in ac_cpp because -g, -O, etc. are not valid cpp options.
ac_cpp='$CXXCPP $CPPFLAGS'
@@ -2192,7 +2224,7 @@ ac_link='${CXX-g++} -o conftest${ac_exeext} $CXXFLAGS $CPPFLAGS $LDFLAGS conftes
cross_compiling=$ac_cv_prog_cxx_cross
cat > conftest.$ac_ext <
@@ -2204,7 +2236,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2208: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2240: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t"""yes"" 1>&6
@@ -2222,7 +2254,7 @@ fi
rm -f conftest*
echo $ac_n "checking whether isnan is declared with debug flags""... $ac_c" 1>&6
-echo "configure:2226: checking whether isnan is declared with debug flags" >&5
+echo "configure:2258: checking whether isnan is declared with debug flags" >&5
ac_ext=C
# CXXFLAGS is not in ac_cpp because -g, -O, etc. are not valid cpp options.
ac_cpp='$CXXCPP $CPPFLAGS'
@@ -2233,7 +2265,7 @@ cross_compiling=$ac_cv_prog_cxx_cross
CXXFLAGS=$CXXFLAGSG
deal_II_isnan_flag=""
cat > conftest.$ac_ext <
@@ -2245,7 +2277,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2249: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2281: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t"""yes"" 1>&6
@@ -2261,7 +2293,7 @@ rm -f conftest*
if test "x$deal_II_isnan_flag" = "x" ; then
cat > conftest.$ac_ext <
@@ -2273,7 +2305,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2277: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2309: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t"""yes"" 1>&6
@@ -2293,7 +2325,7 @@ rm -f conftest*
for testflag in -D_ISOC99_SOURCE -D__EXTENSIONS__ ; do
CXXFLAGS="$CXXFLAGSG $testflag"
cat > conftest.$ac_ext <
@@ -2305,7 +2337,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2309: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2341: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
deal_II_isnan_flag="-DHAVE_ISNAN $testflag"
@@ -2318,7 +2350,7 @@ fi
rm -f conftest*
cat > conftest.$ac_ext <
@@ -2330,7 +2362,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2334: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2366: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
deal_II_isnan_flag="-DHAVE_UNDERSCORE_ISNAN $testflag"
@@ -2352,7 +2384,7 @@ rm -f conftest*
fi
echo $ac_n "checking whether isnan is declared with optimized flags""... $ac_c" 1>&6
-echo "configure:2356: checking whether isnan is declared with optimized flags" >&5
+echo "configure:2388: checking whether isnan is declared with optimized flags" >&5
ac_ext=C
# CXXFLAGS is not in ac_cpp because -g, -O, etc. are not valid cpp options.
ac_cpp='$CXXCPP $CPPFLAGS'
@@ -2363,7 +2395,7 @@ cross_compiling=$ac_cv_prog_cxx_cross
CXXFLAGS=$CXXFLAGSO
deal_II_isnan_flag=""
cat > conftest.$ac_ext <
@@ -2375,7 +2407,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2379: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2411: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t"""yes"" 1>&6
@@ -2391,7 +2423,7 @@ rm -f conftest*
if test "x$deal_II_isnan_flag" = "x" ; then
cat > conftest.$ac_ext <
@@ -2403,7 +2435,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2407: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2439: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t"""yes"" 1>&6
@@ -2423,7 +2455,7 @@ rm -f conftest*
for testflag in -D_ISOC99_SOURCE -D__EXTENSIONS__ ; do
CXXFLAGS="$CXXFLAGSO $testflag"
cat > conftest.$ac_ext <
@@ -2435,7 +2467,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2439: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2471: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
deal_II_isnan_flag="-DHAVE_ISNAN $testflag"
@@ -2448,7 +2480,7 @@ fi
rm -f conftest*
cat > conftest.$ac_ext <
@@ -2460,7 +2492,7 @@ int main() {
; return 0; }
EOF
-if { (eval echo configure:2464: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2496: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
deal_II_isnan_flag="-DHAVE_UNDERSCORE_ISNAN $testflag"
@@ -2495,16 +2527,16 @@ cross_compiling=$ac_cv_prog_cxx_cross
CXXFLAGS="$CXXFLAGSG"
echo $ac_n "checking for consistency of CXXFLAGSG flags""... $ac_c" 1>&6
-echo "configure:2499: checking for consistency of CXXFLAGSG flags" >&5
+echo "configure:2531: checking for consistency of CXXFLAGSG flags" >&5
cat > conftest.$ac_ext <&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2540: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t""yes" 1>&6
@@ -2522,16 +2554,16 @@ rm -f conftest*
CXXFLAGS="$CXXFLAGSO"
echo $ac_n "checking for consistency of CXXFLAGSO flags""... $ac_c" 1>&6
-echo "configure:2526: checking for consistency of CXXFLAGSO flags" >&5
+echo "configure:2558: checking for consistency of CXXFLAGSO flags" >&5
cat > conftest.$ac_ext <&5; (eval $ac_compile) 2>&5; }; then
+if { (eval echo configure:2567: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then
rm -rf conftest*
echo "$ac_t""yes" 1>&6
@@ -2691,6 +2723,7 @@ s%@F77LDFLAGS@%$F77LDFLAGS%g
s%@PERL@%$PERL%g
s%@enablemultithreading@%$enablemultithreading%g
s%@withmultithreading@%$withmultithreading%g
+s%@USE_CONTRIB_HSL@%$USE_CONTRIB_HSL%g
s%@enablemultigrid@%$enablemultigrid%g
s%@kdocdir@%$kdocdir%g
s%@kdocversion@%$kdocversion%g
diff --git a/deal.II/configure.in b/deal.II/configure.in
index c1bfe5ede2..1e0867a60f 100644
--- a/deal.II/configure.in
+++ b/deal.II/configure.in
@@ -545,6 +545,38 @@ if test "$withmultithreading" != no ; then
fi
AC_SUBST(withmultithreading)
+
+
+dnl ------------------------------------------------------------
+dnl Check whether some of the HSL functions have been dropped
+dnl into their respective place in the contrib subdir.
+dnl Check for the following functions to be there:
+dnl MA27: needs files ma27.f
+dnl MA47: needs files ma47.f ma47.dep
+dnl ------------------------------------------------------------
+AC_MSG_CHECKING(for HSL subroutines)
+hsl_subroutines=""
+if test -r contrib/hsl/source/ma27.f ; then
+ hsl_subroutines="$hsl_subroutines MA27"
+ AC_DEFINE(HAVE_HSL_MA27)
+fi
+
+if test -r contrib/hsl/source/ma47.f -a -r contrib/hsl/source/ma47dep.f ; then
+ hsl_subroutines="$hsl_subroutines MA47"
+ AC_DEFINE(HAVE_HSL_MA47)
+fi
+
+AC_MSG_RESULT($hsl_subroutines)
+
+if test "x$hsl_subroutines" != "x" ; then
+ USE_CONTRIB_HSL=yes
+else
+ USE_CONTRIB_HSL=no
+fi
+AC_SUBST(USE_CONTRIB_HSL)
+
+
+
dnl ------------------------------------------------------------
dnl Check for LAPACK and BLAS
dnl Consider this part experimental!
diff --git a/deal.II/contrib/Makefile.in b/deal.II/contrib/Makefile.in
index b70d5ab319..9fa5eb3e2b 100644
--- a/deal.II/contrib/Makefile.in
+++ b/deal.II/contrib/Makefile.in
@@ -2,8 +2,8 @@
# Makefile for the /contrib subdirectory
-# by default do nothing
-default:
+# by default make sublibs
+default: hsl
# targets for the kdoc package
@@ -13,7 +13,13 @@ kdoc-install:
kdoc-clean:
cd kdoc ; $(MAKE) clean
+hsl:
+ cd hsl ; $(MAKE)
# generic targets
clean: kdoc-clean
+ cd hsl ; $(MAKE) clean
+
+
+.PHONY: default kdoc-install kdoc-clean hsl clean
diff --git a/deal.II/contrib/configure.in b/deal.II/contrib/configure.in
index f1f8759b32..4e93042c34 100644
--- a/deal.II/contrib/configure.in
+++ b/deal.II/contrib/configure.in
@@ -27,6 +27,9 @@ cd kdoc/src
cd $DEAL2_CONTRIB_DIR
+dnl Configure the HSL subdir
+AC_CONIG_SUBDIRS(hsl)
+
dnl -------------------------------------------------------------
dnl Output results
diff --git a/deal.II/doc/development/Makefile b/deal.II/doc/development/Makefile
index 30f9c89ccb..8002c8a758 100644
--- a/deal.II/doc/development/Makefile
+++ b/deal.II/doc/development/Makefile
@@ -19,9 +19,11 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
@echo 'PERL=$(PERL)
' >> $@
@echo 'DEAL_II_MAJOR=$(DEAL_II_MAJOR)
' >> $@
@echo 'DEAL_II_MINOR=$(DEAL_II_MINOR)
' >> $@
+ @echo 'USE_CONTRIB_HSL=$(USE_CONTRIB_HSL)
' >> $@
@echo 'lib-path-base=$(lib-path-base)
' >> $@
@echo 'lib-path-lac=$(lib-path-lac)
' >> $@
@echo 'lib-path-deal2=$(lib-path-deal2)
' >> $@
+ @echo 'lib-path-contrib-hsl=$(lib-path-contrib-hsl)
' >> $@
@echo 'LIBPATH=$(LIBPATH)
' >> $@
@echo 'lib-base.o=$(lib-base.o)
' >> $@
@echo 'lib-base.g=$(lib-base.g)
' >> $@
@@ -33,9 +35,11 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
@echo 'lib-deal2-2d.g=$(lib-deal2-2d.g)
' >> $@
@echo 'lib-deal2-3d.o=$(lib-deal2-3d.o)
' >> $@
@echo 'lib-deal2-3d.g=$(lib-deal2-3d.g)
' >> $@
+ @echo 'lib-contrib-hsl=$(lib-contrib-hsl)
' >> $@
@echo 'include-path-base=$(include-path-base)
' >> $@
@echo 'include-path-lac=$(include-path-lac)
' >> $@
@echo 'include-path-deal2=$(include-path-deal2)
' >> $@
+ @echo 'include-path-contrib-hsl=$(include-path-contrib-hsl)
' >> $@
@echo 'INCLUDE=$(INCLUDE)
' >> $@
@echo 'CXXFLAGS.g=$(CXXFLAGS.g)
' >> $@
@echo 'CXXFLAGS.o=$(CXXFLAGS.o)
' >> $@
@@ -43,7 +47,6 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
@echo 'F77FLAGS.o=$(F77FLAGS.o)
' >> $@
@echo 'LDFLAGS=$(LDFLAGS)
' >> $@
@echo 'F77LDFLAGS=$(F77LDFLAGS)
' >> $@
- @echo 'LIBS=$(LIBS)
' >> $@
@echo 'ACE_ROOT=$(ACE_ROOT)
' >> $@
@echo 'lib-ACE=$(lib-ACE)
' >> $@
@echo 'with-multithreading=$(with-multithreading)
' >> $@
diff --git a/deal.II/doc/development/makefiles.1.html b/deal.II/doc/development/makefiles.1.html
index a8a5be3c26..dd09c0351f 100644
--- a/deal.II/doc/development/makefiles.1.html
+++ b/deal.II/doc/development/makefiles.1.html
@@ -85,7 +85,8 @@
Name and version of the Fortran 77 compiler. Possible
names presently include egcs1.1
,
gcc2.95
, gcc2.96
,
- gcc2.97
, SunF77
;
+ gcc2.97
, SunF77
,
+ AIXF77
, MIPSproF77
;
there may be others in the future, if we have access to
other systems
@@ -119,6 +120,15 @@
in programs through a preprocessor variable
+
+
+ USE_CONTRIB_HSL
+
+ yes
indicates whether some subroutines from
+ the HSL have been explicitely installed and shall be used
+ by the library. no
otherwise.
+
+
@@ -147,6 +157,15 @@
+ lib-path-contrib-hsl
+
+ Path to the place where the HSL subroutine functions are
+ linked into a library. Note that this library is only
+ available if it is installed explicitely by the user.
+
+
+
+
LIBPATH
All the library pathes prefixed by -L, i.e. the
@@ -193,6 +212,13 @@
Same for the deal.II libraries in the various dimensions
+
+ lib-contrib-hsl
+
+ Path and name of the library that contains the HSL
+ subroutines explicitely installed.
+
+
@@ -221,6 +247,15 @@
+ include-path-contrib-hsl
+
+ If some HSL files are
+ installed, this path points to includes for the
+ declaration of C bindings to these functions.
+
+
+
+
include-path-deal2
Same for the /deal.II library. Include files are in
@@ -285,15 +320,6 @@
libraries as libg2c
for g77, or vendor libraries
-
-
- LIBS
-
- Additional libraries found by configure and needed to be
- linked with based on the flags given to
- ./configure
, e.g. BLAS or LAPACK
-
-
diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index bee329ab95..4efeb7ce3f 100644
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -26,6 +26,14 @@ documentation, etc.
General
+ -
+ New: There is now some support to include and use routines from the
+ Harwell Subroutine Library.
+
+ (WB 2001/01/30)
+
+
-
New: The ./configure
script now checks for the
existence of a Fortran 77 compiler, and sets its path, some
@@ -68,6 +76,15 @@ documentation, etc.
lac
+ -
+ New: There is now some support to include and use routines from the
+ Harwell Subroutine Library, and support
+ classes for the sparse direct solvers MA27 and MA47.
+
+ (WB 2001/01/30)
+
+
diff --git a/deal.II/lac/Makefile b/deal.II/lac/Makefile
index 71d1640706..501f45a9d5 100644
--- a/deal.II/lac/Makefile
+++ b/deal.II/lac/Makefile
@@ -50,20 +50,20 @@ endif
lib/liblac.a: $(forward-declarations) $(o-files)
- @echo =====lac==============optimized=$(MT)== Linking library: $@
+ @echo =====lac==============optimized==$(MT)== Linking library: $@
@ar ruv $@ $(o-files)
lib/liblac.g.a: $(forward-declarations) $(go-files)
- @echo =====lac==============debug=====$(MT)== Linking library: $@
+ @echo =====lac==============debug======$(MT)== Linking library: $@
@ar ruv $@ $(go-files)
lib/liblac.so: $(forward-declarations) $(o-files)
- @echo =====lac==============optimized=$(MT)== Linking library: $@
+ @echo =====lac==============optimized==$(MT)== Linking library: $@
@$(CC) $(CXXFLAGS.o) -shared -o $@ $(o-files)
lib/liblac.g.so: $(forward-declarations) $(go-files)
- @echo =====lac==============debug=====$(MT)== Linking library: $@
+ @echo =====lac==============debug======$(MT)== Linking library: $@
@$(CC) $(CXXFLAGS.g) -shared -o $@ $(go-files)
diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h
new file mode 100644
index 0000000000..1330e7ab1f
--- /dev/null
+++ b/deal.II/lac/include/lac/sparse_direct.h
@@ -0,0 +1,627 @@
+//---------------------------- sparse_direct.h ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_direct.h ---------------------------
+#ifndef __deal2__sparse_direct_h
+#define __deal2__sparse_direct_h
+
+
+
+#include
+#include
+#include
+
+
+
+/**
+ * This class provides an interface to the sparse direct solver MA27
+ * from the Harwell Subroutine Library. MA27 is a direct solver
+ * specialized for sparse indefinite systems of linear equations and
+ * uses a modified form of Gauss elimination. It is included in the
+ * Harwell Subroutine Library (see
+ * @url{http://www.cse.clrc.ac.uk/Activity/HSL}) and is written in
+ * Fortran. The present class only transforms the data stored in
+ * @ref{SparseMatrix} objects into the form which is required by the
+ * functions resembling MA27, calls these Fortran functions, and
+ * interprets some of the returned values indicating error codes,
+ * etc. It also manages allocation of the right amount of temporary
+ * storage required by these functions.
+ *
+ *
+ * @section3{Interface and Method}
+ *
+ * For the meaning of the three functions @p{initialize},
+ * @p{factorize}, and @p{solve}, as well as for the method used in
+ * MA27, please see the documentation of these functions, which can be
+ * obtained from @url{http://www.cse.clrc.ac.uk/Activity/HSL}. In
+ * practice, one will most often call the second @p{solve} function,
+ * which solves the linear system for a given right hand sidem but one
+ * can as well call the three functions separately if, for example,
+ * one would like to solve the same matrix for several right hand side
+ * vectors; the MA27 solver can do this efficiently, as it computes a
+ * decomposition of the matrix, so that subsequent solves only amount
+ * to a forward-backward substitution which is significantly less
+ * costly than the decomposition process.
+ *
+ *
+ * @section3{Parameters to the constructor}
+ *
+ * The constructor of this class takes several arguments. The meaning
+ * is the following: the MA27 functions require the user to allocate
+ * and pass a certain amount of memory for temporary variables or for
+ * data to be passed to subsequent functions. The sizes of these
+ * arrays are denoted by the variables @p{LIW1}, @p{LIW2}, and @p{LA},
+ * where @p{LIW1} denotes the size of the @p{IW} array in the call to
+ * @p{MA27A}, while @p{LIW2} is the array size in the call to
+ * @p{MA27B}. The documentation of the MA27 functions gives ways to
+ * obtain estimates for their values, e.g. by evaluating values
+ * returned by functions called before. However, the documentation
+ * only states that the values have to be @it{at least as large} as
+ * the estimates, a hint that is not very useful oftentimes (in my
+ * humble opinion, the lack of dynamic memory allocation mechanism is
+ * a good reason not to program in Fortran 77...).
+ *
+ * In our experience, it is often necessary to go beyond the proposed
+ * values (most often for @p{LA}, but also for @p{LIW1}). The first
+ * three parameters of the constructor denote by which factor the
+ * initial estimates shall be increased. The default values are 1.2
+ * (the documentation recommends this value, 1, and 1.5, values which
+ * have often worked for us. Note that the value of @p{LIW} is only
+ * changed in the second call if the recommended value times
+ * @p{LIW_factor_2} is larger than the array size already is from the
+ * call to @p{MA27A}; otherwise, @p{LIW_factor_2} is ignored.
+ *
+ * If the values thus constructed fail to work, we try to restart the
+ * called function with larger values until the calls succeed. The
+ * second triple of values passed to the constructor denotes by which
+ * factor we shall increase the array sizes. If the increment factors
+ * are less than or equal to one, then we only try to call the
+ * respective calls to the functions once and abort by throwing an
+ * error. Note that the @p{MA27C} function writes out an error message
+ * if the value of @p{LA} is too small and gives an indication to
+ * which size it should be increased. However, most often the
+ * indicated value is far too small and can not be relied upon.
+ *
+ * @author Wolfgang Bangerth, 2000, 2001
+ */
+class SparseDirectMA27
+{
+ public:
+ /**
+ * Constructor. See the
+ * documentation of this class
+ * for the meaning of the
+ * parameters to this function.
+ */
+ SparseDirectMA27 (const double LIW_factor_1 = 1.2,
+ const double LIW_factor_2 = 1,
+ const double LA_factor = 1.5,
+ const double LIW_increase_factor_1 = 1.2,
+ const double LIW_increase_factor_2 = 1.2,
+ const double LA_increase_factor = 1.2,
+ const bool suppress_output = true);
+
+ /**
+ * Initialize some data
+ * structures. This function
+ * computes symbolically some
+ * information based on the
+ * sparsity pattern, but does not
+ * actually use the values of the
+ * matrix, so only the sparsity
+ * pattern has to be passed as
+ * argument.
+ */
+ void initialize (const SparsityPattern &sparsity_pattern);
+
+ /**
+ * Actually factorize the
+ * matrix. This function may be
+ * called multiple times for
+ * different matrices, after the
+ * object of this class has been
+ * initialized for a certain
+ * sparsity pattern. You may
+ * therefore save some computing
+ * time if you want to invert
+ * several matrices with the same
+ * sparsity pattern. However,
+ * note that the bulk of the
+ * computing time is actually
+ * spent in the factorization, so
+ * this functionality may not
+ * always be of large benefit.
+ *
+ * If the initialization step has
+ * not been performed yet, then
+ * the @p{initialize} function is
+ * called at the beginning of
+ * this function.
+ */
+ void factorize (const SparseMatrix &matrix);
+
+ /**
+ * Solve for a certain right hand
+ * side vector. This function may
+ * be called multiple times for
+ * different right hand side
+ * vectors after the matrix has
+ * been factorized. This yields a
+ * big saving in computing time,
+ * since the actual solution is
+ * fast, compared to the
+ * factorization of the matrix.
+ *
+ * The solution will be returned
+ * in place of the right hand
+ * side vector.
+ *
+ * If the factorization has not
+ * happened before, strange
+ * things will happen. Note that
+ * we can't actually call the
+ * @p{factorize} function from
+ * here if it has not yet been
+ * called, since we have no
+ * access to the actual matrix.
+ */
+ void solve (Vector &rhs_and_solution) const;
+
+ /**
+ * Call the three functions above
+ * in that order, i.e. perform
+ * the whole solution process for
+ * the given right hand side
+ * vector.
+ *
+ * The solution will be returned
+ * in place of the right hand
+ * side vector.
+ */
+ void solve (const SparseMatrix &matrix,
+ Vector &rhs_and_solution);
+
+ /**
+ * Return an estimate of the
+ * memory used by this class.
+ */
+ unsigned int memory_consumption () const;
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA27AFailed,
+ int,
+ << "The function MA27A failed with an exit code of " << arg1);
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA27BFailed,
+ int,
+ << "The function MA27B failed with an exit code of " << arg1);
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA27CFailed,
+ int,
+ << "The function MA27C failed with an exit code of " << arg1);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInitializeAlreadyCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcFactorizeNotCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcDifferentSparsityPatterns);
+
+ private:
+ /**
+ * Store the three values passed
+ * to the cinstructor. See the
+ * documentation of this class
+ * for the meaning of these
+ * variables.
+ */
+ const double LIW_factor_1;
+ const double LIW_factor_2;
+ const double LA_factor;
+
+ /**
+ * Increase factors in case a
+ * call to a function fails.
+ */
+ const double LIW_increase_factor_1;
+ const double LIW_increase_factor_2;
+ const double LA_increase_factor;
+
+ /**
+ * Flags storing whether the
+ * first two functions have
+ * already been called.
+ */
+ bool initialize_called;
+ bool factorize_called;
+
+ /**
+ * Store a pointer to the
+ * sparsity pattern, to make sure
+ * that we use the same thing for
+ * all calls.
+ */
+ SmartPointer sparsity_pattern;
+
+ /**
+ * Number of nonzero elements in
+ * the sparsity pattern on and
+ * above the diagonal.
+ */
+ unsigned int n_nonzero_elements;
+
+ /**
+ * Arrays holding row and column
+ * indices.
+ */
+ vector row_numbers;
+ vector column_numbers;
+
+ /**
+ * Array to hold the matrix
+ * elements, and later the
+ * elements of the factors.
+ */
+ vector A;
+
+ /**
+ * Length of the @p{A} array.
+ */
+ unsigned int LA;
+
+ /**
+ * Scratch arrays and variables
+ * used by the MA27 functions. We
+ * keep to the names introduced
+ * in the documentation of these
+ * functions, in all uppercase
+ * letters as is usual in
+ * Fortran.
+ */
+ unsigned int LIW;
+ vector IW;
+ vector IKEEP;
+ vector IW1;
+
+ unsigned int NSTEPS;
+ unsigned int MAXFRT;
+
+ /**
+ * Two values that live inside a
+ * COMMON block of the Fortran
+ * code and are mirrored at these
+ * locations. They are used to
+ * transport information about
+ * the required length of arrays
+ * from the Fortran functions to
+ * the outside world.
+ */
+ unsigned int NRLNEC;
+ unsigned int NIRNEC;
+
+ /**
+ * Flag indicating the level of
+ * output desired and returning
+ * error values if error occured.
+ */
+ int IFLAG;
+
+ /**
+ * Fill the @p{A} array from the
+ * symmetric part of the given
+ * matrix.
+ */
+ void fill_A (const SparseMatrix &matrix);
+};
+
+
+
+
+
+
+/**
+ * This class provides an interface to the sparse direct solver MA47
+ * from the Harwell Subroutine Library. MA47 is a direct solver
+ * specialized for sparse indefinite systems of linear equations and
+ * uses a frontal elimination method. It is included in the Harwell
+ * Subroutine Library (see
+ * @url{http://www.cse.clrc.ac.uk/Activity/HSL}) and is written in
+ * Fortran. The present class only transforms the data stored in
+ * @ref{SparseMatrix} objects into the form which is required by the
+ * functions resembling MA47, calls these Fortran functions, and
+ * interprets some of the returned values indicating error codes,
+ * etc. It also manages allocation of the right amount of temporary
+ * storage required by these functions.
+ *
+ *
+ * @section3{Interface and Method}
+ *
+ * For the meaning of the three functions @p{initialize},
+ * @p{factorize}, and @p{solve}, as well as for the method used in
+ * MA47, please see the documentation of these functions, which can be
+ * obtained from @url{http://www.cse.clrc.ac.uk/Activity/HSL}. In
+ * practice, one will most often call the second @p{solve} function,
+ * which solves the linear system for a given right hand sidem but one
+ * can as well call the three functions separately if, for example,
+ * one would like to solve the same matrix for several right hand side
+ * vectors; the MA47 solver can do this efficiently, as it computes a
+ * decomposition of the matrix, so that subsequent solves only amount
+ * to a forward-backward substitution which is significantly less
+ * costly than the decomposition process.
+ *
+ *
+ * @section3{Parameters to the constructor}
+ *
+ * The constructor of this class takes several arguments. Their
+ * meaning is equivalent to those of the constructor of the
+ * @ref{SparseDirectMA27} class; see there for more information.
+ *
+ * @author Wolfgang Bangerth, 2000, 2001
+ */
+class SparseDirectMA47
+{
+ public:
+ /**
+ * Constructor. See the
+ * documentation of this class
+ * for the meaning of the
+ * parameters to this function.
+ *
+ * This function already calls
+ * the initialization function
+ * @p{MA47ID} to set up some
+ * values.
+ */
+ SparseDirectMA47 (const double LIW_factor_1 = 1.4,
+ const double LIW_factor_2 = 1,
+ const double LA_factor = 3,
+ const double LIW_increase_factor_1 = 1.2,
+ const double LIW_increase_factor_2 = 1.2,
+ const double LA_increase_factor = 1.2,
+ const bool suppress_output = true);
+
+ /**
+ * Initialize some data
+ * structures. This function
+ * computes symbolically some
+ * information based on the
+ * sparsity pattern, but does not
+ * actually use the values of the
+ * matrix, so only the sparsity
+ * pattern has to be passed as
+ * argument.
+ *
+ * Since the MA47 solver requires
+ * us to omit zero-entries in the
+ * matrix (even if they are in
+ * the sparsity pattern), we have
+ * to actually use the matrix
+ * here, as opposed to the MA27
+ * solver that only required the
+ * sparsity pattern.
+ */
+ void initialize (const SparseMatrix &matrix);
+
+ /**
+ * Actually factorize the
+ * matrix. Unlike for the MA27
+ * solver, this function may not
+ * be called multiple times for
+ * different matrices, since we
+ * have eliminated entries from
+ * the sparsity pattern where
+ * matrix entries happen to be
+ * zero. Since this is likely to
+ * change between matrices
+ * although they have the same
+ * sparsity pattern.
+ *
+ * If the initialization step has
+ * not been performed yet, then
+ * the @p{initialize} function is
+ * called at the beginning of
+ * this function.
+ */
+ void factorize (const SparseMatrix &matrix);
+
+ /**
+ * Solve for a certain right hand
+ * side vector. This function may
+ * be called multiple times for
+ * different right hand side
+ * vectors after the matrix has
+ * been factorized. This yields a
+ * big saving in computing time,
+ * since the actual solution is
+ * fast, compared to the
+ * factorization of the matrix.
+ *
+ * The solution will be returned
+ * in place of the right hand
+ * side vector.
+ *
+ * If the factorization has not
+ * happened before, strange
+ * things will happen. Note that
+ * we can't actually call the
+ * @p{factorize} function from
+ * here if it has not yet been
+ * called, since we have no
+ * access to the actual matrix.
+ */
+ void solve (Vector &rhs_and_solution);
+
+ /**
+ * Call the three functions above
+ * in that order, i.e. perform
+ * the whole solution process for
+ * the given right hand side
+ * vector.
+ *
+ * The solution will be returned
+ * in place of the right hand
+ * side vector.
+ */
+ void solve (const SparseMatrix &matrix,
+ Vector &rhs_and_solution);
+
+ /**
+ * Return an estimate of the
+ * memory used by this class.
+ */
+ unsigned int memory_consumption () const;
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA47AFailed,
+ int,
+ << "The function MA47A failed with an exit code of " << arg1);
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA47BFailed,
+ int,
+ << "The function MA47B failed with an exit code of " << arg1);
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcMA47CFailed,
+ int,
+ << "The function MA47C failed with an exit code of " << arg1);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInitializeAlreadyCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcFactorizeNotCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcCantFactorizeAgain);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcDifferentMatrices);
+
+ private:
+ /**
+ * Store the three values passed
+ * to the cinstructor. See the
+ * documentation of this class
+ * for the meaning of these
+ * variables.
+ */
+ const double LIW_factor_1;
+ const double LIW_factor_2;
+ const double LA_factor;
+
+ /**
+ * Increase factors in case a
+ * call to a function fails.
+ */
+ const double LIW_increase_factor_1;
+ const double LIW_increase_factor_2;
+ const double LA_increase_factor;
+
+ /**
+ * Flags storing whether the
+ * first two functions have
+ * already been called.
+ */
+ bool initialize_called;
+ bool factorize_called;
+
+ /**
+ * Store a pointer to the matrix,
+ * to make sure that we use the
+ * same thing for all calls.
+ */
+ SmartPointer > matrix;
+
+ /**
+ * Number of nonzero elements in
+ * the sparsity pattern on and
+ * above the diagonal.
+ */
+ unsigned int n_nonzero_elements;
+
+ /**
+ * Control values set by @p{MA47ID}.
+ */
+ double CNTL[2];
+ unsigned int ICNTL[7];
+
+ /**
+ * Info field filled by the MA47
+ * functions and (partially) used
+ * for subsequent MA47 calls.
+ */
+ int INFO[24];
+
+ /**
+ * Arrays holding row and column
+ * indices.
+ */
+ vector row_numbers;
+ vector column_numbers;
+
+ /**
+ * Array to hold the matrix
+ * elements, and later the
+ * elements of the factors.
+ */
+ vector A;
+
+ /**
+ * Length of the @p{A} array.
+ */
+ unsigned int LA;
+
+ /**
+ * Scratch arrays and variables
+ * used by the MA47 functions. We
+ * keep to the names introduced
+ * in the documentation of these
+ * functions, in all uppercase
+ * letters as is usual in
+ * Fortran.
+ */
+ unsigned int LIW;
+ vector IW;
+ vector KEEP;
+ vector IW1;
+
+ /**
+ * Fill the @p{A} array from the
+ * symmetric part of the given
+ * matrix.
+ */
+ void fill_A (const SparseMatrix &matrix);
+};
+
+
+
+
+#endif
diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc
new file mode 100644
index 0000000000..b0879d5c19
--- /dev/null
+++ b/deal.II/lac/source/sparse_direct.cc
@@ -0,0 +1,855 @@
+//---------------------------- sparse_direct.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_direct.cc ---------------------------
+
+
+#include
+#include
+#include
+#include
+#include
+
+
+// if the HSL functions are not there, define them empty and throw an
+// exception
+#ifndef HAVE_HSL_MA27
+namespace HSL
+{
+ namespace MA27
+ {
+ extern "C"
+ void ma27ad_ (const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ unsigned int *,
+ const unsigned int *,
+ unsigned int *,
+ unsigned int *,
+ unsigned int *,
+ int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C"
+ void ma27bd_ (const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ double *,
+ const unsigned int *,
+ unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ unsigned int *,
+ unsigned int *,
+ int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C"
+ void ma27cd_ (const unsigned int *,
+ const double *,
+ const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ double *,
+ const unsigned int *,
+ double *,
+ const unsigned int *,
+ const unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C" void ma27x1_ (unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C" void ma27x2_ (unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C" void ma27x3_ (const unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+ };
+};
+#endif // ifndef HAVE_HSL_MA27
+
+
+#ifndef HAVE_HSL_MA47
+namespace HSL
+{
+ namespace MA47
+ {
+ extern "C"
+ void ma47id_ (double *,
+ unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C"
+ void ma47ad_ (const unsigned int *,
+ const unsigned int *,
+ unsigned int *,
+ unsigned int *,
+ unsigned int *,
+ const unsigned int *,
+ unsigned int *,
+ const unsigned int *,
+ double *,
+ int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C"
+ void ma47bd_ (const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ double *,
+ const unsigned int *,
+ unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ const double *,
+ const unsigned int *,
+ unsigned int *,
+ double *,
+ int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+
+
+ extern "C"
+ void ma47cd_ (const unsigned int *,
+ const double *,
+ const unsigned int *,
+ const unsigned int *,
+ const unsigned int *,
+ double *,
+ double *,
+ unsigned int *,
+ const unsigned int *)
+ {
+ AssertThrow (false, ExcMessage("You can only use the HSL functions after putting "
+ "the respective files in the right place, "
+ "re-configuring the library and re-building it!"));
+ };
+ };
+};
+#endif // ifndef HAVE_HSL_MA47
+
+
+
+
+
+/* -------------------------- MA27 ---------------------------- */
+
+
+SparseDirectMA27::SparseDirectMA27 (const double LIW_factor_1,
+ const double LIW_factor_2,
+ const double LA_factor,
+ const double LIW_increase_factor_1,
+ const double LIW_increase_factor_2,
+ const double LA_increase_factor,
+ const bool suppress_output) :
+ LIW_factor_1 (LIW_factor_1),
+ LIW_factor_2 (LIW_factor_2),
+ LA_factor (LA_factor),
+ LIW_increase_factor_1 (LIW_increase_factor_1),
+ LIW_increase_factor_2 (LIW_increase_factor_2),
+ LA_increase_factor (LA_increase_factor),
+ initialize_called (false),
+ factorize_called (false),
+ sparsity_pattern (0)
+{
+ // suppress error output if
+ // requested
+ if (suppress_output)
+ {
+ const unsigned int LP = 0;
+ HSL::MA27::ma27x3_ (&LP);
+ };
+};
+
+
+
+void
+SparseDirectMA27::initialize (const SparsityPattern &sp)
+{
+ Assert (initialize_called == false, ExcInitializeAlreadyCalled());
+
+ sparsity_pattern = &sp;
+
+ const unsigned int n_rows = sparsity_pattern->n_rows();
+ const unsigned int *rowstart_indices = sparsity_pattern->get_rowstart_indices();
+ const unsigned int *col_nums = sparsity_pattern->get_column_numbers();
+
+ // first count number of nonzero
+ // elements in the upper right
+ // part. the matrix is symmetric,
+ // so this suffices
+ n_nonzero_elements = 0;
+ for (unsigned int row=0; row((2*n_nonzero_elements + 3*n_rows + 1) *
+ LIW_factor_1);
+ IW.resize (LIW);
+ IKEEP.resize (3*n_rows);
+ IW1.resize (2*n_rows);
+
+ // no output please
+ IFLAG = 0;
+
+ // loop until memory requirements
+ // are satisfied or we are not
+ // allowed to allocate more memory
+ // no more
+ bool call_succeeded = true;
+ do
+ {
+ HSL::MA27::ma27ad_(&n_rows, &n_nonzero_elements,
+ &row_numbers[0], &column_numbers[0],
+ &IW[0], &LIW, &IKEEP[0],
+ &IW1[0], &NSTEPS, &IFLAG);
+ call_succeeded = (IFLAG==0);
+
+ // if enough memory or no
+ // increase allowed: exit loop
+ if (call_succeeded || (LIW_increase_factor_1 <= 1))
+ break;
+
+ // otherwise: increase LIW and retry
+ LIW = static_cast(LIW * LIW_increase_factor_1);
+ IW.resize (LIW);
+ }
+ while (true);
+
+ // if we were not allowed to
+ // allocate more memory, then throw
+ // an exception
+ if (!call_succeeded)
+ throw ExcMA27AFailed(IFLAG);
+
+ // catch returned values from the
+ // COMMON block. we need these
+ // values in order to set array
+ // sizes in the next function
+ HSL::MA27::ma27x1_(&NRLNEC);
+ HSL::MA27::ma27x2_(&NIRNEC);
+
+ // note that we have already been
+ // in this function
+ initialize_called = true;
+};
+
+
+
+void
+SparseDirectMA27::factorize (const SparseMatrix &matrix)
+{
+ // if necessary, initialize process
+ if (initialize_called == false)
+ initialize (matrix.get_sparsity_pattern());
+
+ // make sure the sparsity patterns
+ // are the same
+ Assert (sparsity_pattern == &matrix.get_sparsity_pattern(),
+ ExcDifferentSparsityPatterns());
+
+
+ // set LA and fill the A array of
+ // values
+ LA=static_cast(NRLNEC * LA_factor);
+ A.resize (LA);
+ fill_A (matrix);
+
+ // if necessary extend IW
+ if (LIW < NIRNEC * LIW_factor_2)
+ {
+ LIW = static_cast(NIRNEC * LIW_factor_2);
+ IW.resize (LIW);
+ };
+
+ const unsigned int n_rows = matrix.get_sparsity_pattern().n_rows();
+
+ // loop until memory requirements
+ // are satisfied or we are not
+ // allowed to allocate more memory
+ // no more
+ bool call_succeeded = true;
+ do
+ {
+ HSL::MA27::ma27bd_(&n_rows, &n_nonzero_elements,
+ &row_numbers[0], &column_numbers[0],
+ &A[0], &LA,
+ &IW[0], &LIW, &IKEEP[0], &NSTEPS, &MAXFRT,
+ &IW1[0], &IFLAG);
+ call_succeeded = (IFLAG==0);
+
+ // if enough memory or no
+ // increase allowed: exit loop
+ if (call_succeeded)
+ break;
+
+ // otherwise: increase LIW or
+ // LA if that is allowed and
+ // retry
+ switch (IFLAG)
+ {
+ case -3:
+ {
+ if (LIW_increase_factor_2 <= 1)
+ goto exit_loop;
+
+ LIW = static_cast(LIW * LIW_increase_factor_2);
+ IW.resize (LIW);
+ break;
+ };
+
+ case -4:
+ {
+ if (LA_increase_factor <= 1)
+ goto exit_loop;
+ // increase A. note that
+ // since the function has
+ // already part of the
+ // array @p{A}, we have
+ // to re-fill it with the
+ // original values. minor
+ // clue: since the old
+ // entries are no more
+ // needed, we can discard
+ // them; we use this to
+ // first release all
+ // memory (through the
+ // call to @p{swap} and
+ // the subsequent call to
+ // the destructor of the
+ // @p{tmp} object) and
+ // only then re-allocate
+ // it. If we called
+ // @p{resize} directly,
+ // this would first
+ // allocate more memory,
+ // then copy the old
+ // contents, and only
+ // then release the old
+ // memory, but keeping
+ // both memory regions at
+ // the same time could
+ // sometimes be more than
+ // we can do, leading to
+ // an exception on the
+ // allocation.
+ cout << "<*>" << flush;
+
+ LA = static_cast(LA * LA_increase_factor);
+ if (true)
+ {
+ vector tmp;
+ A.swap (tmp);
+ };
+
+ A.resize (LA);
+ fill_A (matrix);
+
+ break;
+ };
+
+ // ups, other return
+ // value, don't know
+ // what to do here
+ default:
+ throw ExcMA27BFailed(IFLAG);
+ };
+ continue;
+
+ exit_loop:
+ break;
+ }
+ while (true);
+
+ if (!call_succeeded)
+ throw ExcMA27BFailed(IFLAG);
+
+ // note that we have been here
+ // already
+ factorize_called = true;
+};
+
+
+
+void
+SparseDirectMA27::solve (Vector &rhs_and_solution) const
+{
+ Assert (factorize_called == true, ExcFactorizeNotCalled());
+
+ const unsigned int n_rows = rhs_and_solution.size();
+ vector W(MAXFRT);
+ HSL::MA27::ma27cd_(&n_rows, &A[0], &LA,
+ &IW[0], &LIW, &W[0], &MAXFRT,
+ &rhs_and_solution(0), &IW1[0], &NSTEPS);
+};
+
+
+
+void
+SparseDirectMA27::solve (const SparseMatrix &matrix,
+ Vector &rhs_and_solution)
+{
+ initialize (matrix.get_sparsity_pattern());
+ factorize (matrix);
+ solve (rhs_and_solution);
+};
+
+
+
+unsigned int
+SparseDirectMA27::memory_consumption () const
+{
+ return (sizeof(*this) +
+ MemoryConsumption::memory_consumption (row_numbers) +
+ MemoryConsumption::memory_consumption (column_numbers) +
+ MemoryConsumption::memory_consumption (A) +
+ MemoryConsumption::memory_consumption (IW) +
+ MemoryConsumption::memory_consumption (IKEEP) +
+ MemoryConsumption::memory_consumption (IW1));
+};
+
+
+
+void
+SparseDirectMA27::fill_A (const SparseMatrix &matrix)
+{
+
+ const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
+
+ const unsigned int n_rows = sparsity_pattern.n_rows();
+ const unsigned int *rowstart_indices = sparsity_pattern.get_rowstart_indices();
+ const unsigned int *col_nums = sparsity_pattern.get_column_numbers();
+
+ unsigned int global_index = 0;
+ for (unsigned int row=0; row &m)
+{
+ Assert (initialize_called == false, ExcInitializeAlreadyCalled());
+
+ matrix = &m;
+ const SparsityPattern &sparsity_pattern = matrix->get_sparsity_pattern();
+
+ const unsigned int n_rows = sparsity_pattern.n_rows();
+ const unsigned int *rowstart_indices = sparsity_pattern.get_rowstart_indices();
+ const unsigned int *col_nums = sparsity_pattern.get_column_numbers();
+
+ // first count number of nonzero
+ // elements in the upper right
+ // part. the matrix is symmetric,
+ // so this suffices
+ n_nonzero_elements = 0;
+ for (unsigned int row=0; row((2*n_nonzero_elements + 5*n_rows + 4) *
+ LIW_factor_1);
+ IW.resize (LIW);
+ KEEP.resize (n_nonzero_elements + 5*n_rows + 2);
+
+ // declare output info fields
+ double RINFO[4];
+ bool call_succeeded;
+ do
+ {
+ HSL::MA47::ma47ad_(&n_rows, &n_nonzero_elements,
+ &row_numbers[0], &column_numbers[0],
+ &IW[0], &LIW, &KEEP[0],
+ &ICNTL[0], &RINFO[0], &INFO[0]);
+ call_succeeded = (INFO[0] == 0);
+
+ // if enough memory or no
+ // increase allowed: exit loop
+ if (call_succeeded || (LIW_increase_factor_1 <= 1))
+ break;
+
+ // otherwise: increase LIW and retry
+ LIW = static_cast(LIW * LIW_increase_factor_1);
+ IW.resize (LIW);
+ }
+ while (true);
+
+ if (!call_succeeded)
+ throw ExcMA47AFailed(INFO[0]);
+
+ // note that we have already been
+ // in this function
+ initialize_called = true;
+};
+
+
+
+void
+SparseDirectMA47::factorize (const SparseMatrix &m)
+{
+ Assert (factorize_called == false, ExcCantFactorizeAgain());
+
+ // if necessary, initialize process
+ if (initialize_called == false)
+ initialize (m);
+
+ // make sure the matrices
+ // are the same
+ Assert (matrix == &m, ExcDifferentMatrices());
+
+
+ // set LA and fill the A array of
+ // values
+ LA=static_cast(INFO[5] * LA_factor);
+ A.resize (LA);
+ fill_A (m);
+
+ // if necessary extend IW
+ if (LIW < INFO[6] * LIW_factor_2)
+ {
+ LIW = static_cast(INFO[6] * LIW_factor_2);
+ IW.resize (LIW);
+ };
+
+ const unsigned int n_rows = m.get_sparsity_pattern().n_rows();
+ IW1.resize (2*n_rows+2);
+
+ // output info flags
+ double RINFO[4];
+
+ bool call_succeeded;
+ do
+ {
+ HSL::MA47::ma47bd_(&n_rows, &n_nonzero_elements, &column_numbers[0],
+ &A[0], &LA,
+ &IW[0], &LIW, &KEEP[0], &CNTL[0], &ICNTL[0],
+ &IW1[0], &RINFO[0], &INFO[0]);
+ call_succeeded = (INFO[0] == 0);
+
+ // if enough memory or no
+ // increase allowed: exit loop
+ if (call_succeeded)
+ break;
+
+ // otherwise: increase LIW or
+ // LA if that is allowed and
+ // retry
+ switch (INFO[0])
+ {
+ case -3:
+ {
+ if (LIW_increase_factor_2 <= 1)
+ goto exit_loop;
+
+ LIW = static_cast(LIW * LIW_increase_factor_2);
+ IW.resize (LIW);
+ break;
+ };
+
+ case -4:
+ {
+ if (LA_increase_factor <= 1)
+ goto exit_loop;
+ // increase A. note that
+ // since the function has
+ // already part of the
+ // array @p{A}, we have
+ // to re-fill it with the
+ // original values. minor
+ // clue: since the old
+ // entries are no more
+ // needed, we can discard
+ // them; we use this to
+ // first release all
+ // memory (through the
+ // call to @p{swap} and
+ // the subsequent call to
+ // the destructor of the
+ // @p{tmp} object) and
+ // only then re-allocate
+ // it. If we called
+ // @p{resize} directly,
+ // this would first
+ // allocate more memory,
+ // then copy the old
+ // contents, and only
+ // then release the old
+ // memory, but keeping
+ // both memory regions at
+ // the same time could
+ // sometimes be more than
+ // we can do, leading to
+ // an exception on the
+ // allocation.
+ cout << "<*>" << flush;
+
+ LA = static_cast(LA * LA_increase_factor);
+ if (true)
+ {
+ vector tmp;
+ A.swap (tmp);
+ };
+
+ A.resize (LA);
+ fill_A (m);
+
+ break;
+ };
+
+ // ups, other return
+ // value, don't know
+ // what to do here
+ default:
+ throw ExcMA47BFailed(INFO[0]);
+ };
+ continue;
+
+ exit_loop:
+ break;
+ }
+ while (true);
+
+ if (!call_succeeded)
+ throw ExcMA47BFailed(INFO[0]);
+
+ // note that we have been here
+ // already
+ factorize_called = true;
+};
+
+
+
+void
+SparseDirectMA47::solve (Vector &rhs_and_solution)
+{
+ Assert (factorize_called == true, ExcFactorizeNotCalled());
+
+ const unsigned int n_rows = rhs_and_solution.size();
+ vector W(n_rows);
+ HSL::MA47::ma47cd_(&n_rows, &A[0], &LA,
+ &IW[0], &LIW, &W[0],
+ &rhs_and_solution(0), &IW1[0], &ICNTL[0]);
+};
+
+
+
+void
+SparseDirectMA47::solve (const SparseMatrix &matrix,
+ Vector &rhs_and_solution)
+{
+ initialize (matrix);
+ factorize (matrix);
+ solve (rhs_and_solution);
+};
+
+
+
+unsigned int
+SparseDirectMA47::memory_consumption () const
+{
+ return (sizeof(*this) +
+ MemoryConsumption::memory_consumption (row_numbers) +
+ MemoryConsumption::memory_consumption (column_numbers) +
+ MemoryConsumption::memory_consumption (A) +
+ MemoryConsumption::memory_consumption (IW) +
+ MemoryConsumption::memory_consumption (KEEP) +
+ MemoryConsumption::memory_consumption (IW1));
+};
+
+
+
+void
+SparseDirectMA47::fill_A (const SparseMatrix &matrix)
+{
+ const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
+
+ const unsigned int n_rows = sparsity_pattern.n_rows();
+ const unsigned int *rowstart_indices = sparsity_pattern.get_rowstart_indices();
+ const unsigned int *col_nums = sparsity_pattern.get_column_numbers();
+
+ unsigned int global_index = 0;
+ for (unsigned int row=0; row