From b9c7852abd1f83263d070564274f001baf4ee5cd Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 30 Jan 2001 16:38:28 +0000 Subject: [PATCH] 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 --- deal.II/Makefile | 14 +- deal.II/common/Make.global_options.in | 46 +- deal.II/configure | 109 ++- deal.II/configure.in | 32 + deal.II/contrib/Makefile.in | 10 +- deal.II/contrib/configure.in | 3 + deal.II/doc/development/Makefile | 5 +- deal.II/doc/development/makefiles.1.html | 46 +- deal.II/doc/news/2001/c-3-1.html | 17 + deal.II/lac/Makefile | 8 +- deal.II/lac/include/lac/sparse_direct.h | 627 +++++++++++++++++ deal.II/lac/source/sparse_direct.cc | 855 +++++++++++++++++++++++ 12 files changed, 1702 insertions(+), 70 deletions(-) create mode 100644 deal.II/lac/include/lac/sparse_direct.h create mode 100644 deal.II/lac/source/sparse_direct.cc 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

      +
    1. + New: There is now some support to include and use routines from the + Harwell Subroutine Library. +
      + (WB 2001/01/30) +

      +
    2. 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

        +
      1. + 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