@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. ="
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
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)
.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
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
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)
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)
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
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
ac_save_LIBS="$LIBS"
LIBS="-lblas $LIBS"
cat > conftest.$ac_ext <<EOF
-#line 1958 "configure"
+#line 1990 "configure"
#include "confdefs.h"
int main() {
main()
; return 0; }
EOF
-if { (eval echo configure:1965: \"$ac_link\") 1>&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
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
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
ac_save_LIBS="$LIBS"
LIBS="-llapack $LIBS"
cat > conftest.$ac_ext <<EOF
-#line 2021 "configure"
+#line 2053 "configure"
#include "confdefs.h"
/* Override any gcc2 internal prototype to avoid an error. */
#ifdef __cplusplus
dgels_()
; return 0; }
EOF
-if { (eval echo configure:2035: \"$ac_link\") 1>&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
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
# 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
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
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'
cross_compiling=$ac_cv_prog_cxx_cross
cat > conftest.$ac_ext <<EOF
-#line 2196 "configure"
+#line 2228 "configure"
#include "confdefs.h"
#include <sys/resource.h>
; 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
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'
CXXFLAGS=$CXXFLAGSG
deal_II_isnan_flag=""
cat > conftest.$ac_ext <<EOF
-#line 2237 "configure"
+#line 2269 "configure"
#include "confdefs.h"
#include <cmath>
; 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
if test "x$deal_II_isnan_flag" = "x" ; then
cat > conftest.$ac_ext <<EOF
-#line 2265 "configure"
+#line 2297 "configure"
#include "confdefs.h"
#include <cmath>
; 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
for testflag in -D_ISOC99_SOURCE -D__EXTENSIONS__ ; do
CXXFLAGS="$CXXFLAGSG $testflag"
cat > conftest.$ac_ext <<EOF
-#line 2297 "configure"
+#line 2329 "configure"
#include "confdefs.h"
#include <cmath>
; 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"
rm -f conftest*
cat > conftest.$ac_ext <<EOF
-#line 2322 "configure"
+#line 2354 "configure"
#include "confdefs.h"
#include <cmath>
; 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"
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'
CXXFLAGS=$CXXFLAGSO
deal_II_isnan_flag=""
cat > conftest.$ac_ext <<EOF
-#line 2367 "configure"
+#line 2399 "configure"
#include "confdefs.h"
#include <cmath>
; 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
if test "x$deal_II_isnan_flag" = "x" ; then
cat > conftest.$ac_ext <<EOF
-#line 2395 "configure"
+#line 2427 "configure"
#include "confdefs.h"
#include <cmath>
; 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
for testflag in -D_ISOC99_SOURCE -D__EXTENSIONS__ ; do
CXXFLAGS="$CXXFLAGSO $testflag"
cat > conftest.$ac_ext <<EOF
-#line 2427 "configure"
+#line 2459 "configure"
#include "confdefs.h"
#include <cmath>
; 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"
rm -f conftest*
cat > conftest.$ac_ext <<EOF
-#line 2452 "configure"
+#line 2484 "configure"
#include "confdefs.h"
#include <cmath>
; 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"
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 <<EOF
-#line 2501 "configure"
+#line 2533 "configure"
#include "confdefs.h"
int main() {
;
; return 0; }
EOF
-if { (eval echo configure:2508: \"$ac_compile\") 1>&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
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 <<EOF
-#line 2528 "configure"
+#line 2560 "configure"
#include "confdefs.h"
int main() {
;
; return 0; }
EOF
-if { (eval echo configure:2535: \"$ac_compile\") 1>&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
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
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!
# Makefile for the /contrib subdirectory
-# by default do nothing
-default:
+# by default make sublibs
+default: hsl
# targets for the kdoc package
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
cd $DEAL2_CONTRIB_DIR
+dnl Configure the HSL subdir
+AC_CONIG_SUBDIRS(hsl)
+
dnl -------------------------------------------------------------
dnl Output results
@echo '<li><code>PERL=$(PERL)</code>' >> $@
@echo '<li><code>DEAL_II_MAJOR=$(DEAL_II_MAJOR)</code>' >> $@
@echo '<li><code>DEAL_II_MINOR=$(DEAL_II_MINOR)</code>' >> $@
+ @echo '<li><code>USE_CONTRIB_HSL=$(USE_CONTRIB_HSL)</code>' >> $@
@echo '<li><code>lib-path-base=$(lib-path-base)</code>' >> $@
@echo '<li><code>lib-path-lac=$(lib-path-lac)</code>' >> $@
@echo '<li><code>lib-path-deal2=$(lib-path-deal2)</code>' >> $@
+ @echo '<li><code>lib-path-contrib-hsl=$(lib-path-contrib-hsl)</code>' >> $@
@echo '<li><code>LIBPATH=$(LIBPATH)</code>' >> $@
@echo '<li><code>lib-base.o=$(lib-base.o)</code>' >> $@
@echo '<li><code>lib-base.g=$(lib-base.g)</code>' >> $@
@echo '<li><code>lib-deal2-2d.g=$(lib-deal2-2d.g)</code>' >> $@
@echo '<li><code>lib-deal2-3d.o=$(lib-deal2-3d.o)</code>' >> $@
@echo '<li><code>lib-deal2-3d.g=$(lib-deal2-3d.g)</code>' >> $@
+ @echo '<li><code>lib-contrib-hsl=$(lib-contrib-hsl)</code>' >> $@
@echo '<li><code>include-path-base=$(include-path-base)</code>' >> $@
@echo '<li><code>include-path-lac=$(include-path-lac)</code>' >> $@
@echo '<li><code>include-path-deal2=$(include-path-deal2)</code>' >> $@
+ @echo '<li><code>include-path-contrib-hsl=$(include-path-contrib-hsl)</code>' >> $@
@echo '<li><code>INCLUDE=$(INCLUDE)</code>' >> $@
@echo '<li><code>CXXFLAGS.g=$(CXXFLAGS.g)</code>' >> $@
@echo '<li><code>CXXFLAGS.o=$(CXXFLAGS.o)</code>' >> $@
@echo '<li><code>F77FLAGS.o=$(F77FLAGS.o)</code>' >> $@
@echo '<li><code>LDFLAGS=$(LDFLAGS)</code>' >> $@
@echo '<li><code>F77LDFLAGS=$(F77LDFLAGS)</code>' >> $@
- @echo '<li><code>LIBS=$(LIBS)</code>' >> $@
@echo '<li><code>ACE_ROOT=$(ACE_ROOT)</code>' >> $@
@echo '<li><code>lib-ACE=$(lib-ACE)</code>' >> $@
@echo '<li><code>with-multithreading=$(with-multithreading)</code>' >> $@
Name and version of the Fortran 77 compiler. Possible
names presently include <code>egcs1.1</code>,
<code>gcc2.95</code>, <code>gcc2.96</code>,
- <code>gcc2.97</code>, <code>SunF77</code>;
+ <code>gcc2.97</code>, <code>SunF77</code>,
+ <code>AIXF77</code>, <code>MIPSproF77</code>;
there may be others in the future, if we have access to
other systems
</p>
in programs through a preprocessor variable
</p>
</dd>
+
+
+ <dt> <code>USE_CONTRIB_HSL</code> </dt>
+ <dd> <p>
+ <code>yes</code> indicates whether some subroutines from
+ the HSL have been explicitely installed and shall be used
+ by the library. <code>no</code> otherwise.
+ </p>
+ </dd>
</dl>
</dd>
+ <dt> <code>lib-path-contrib-hsl</code> </dt>
+ <dd> <p>
+ 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.
+ </p>
+ </dd>
+
+
<dt> <code>LIBPATH</code> </dt>
<dd> <p>
All the library pathes prefixed by -L, i.e. the
Same for the deal.II libraries in the various dimensions
</p>
</dd>
+
+ <dt> <code>lib-contrib-hsl</code> </dt>
+ <dd> <p>
+ Path and name of the library that contains the HSL
+ subroutines explicitely installed.
+ </p>
+ </dd>
</dl>
</dd>
+ <dt> <code>include-path-contrib-hsl</code> </dt>
+ <dd> <p>
+ If some HSL files are
+ installed, this path points to includes for the
+ declaration of C bindings to these functions.
+ </p>
+ </dd>
+
+
<dt> <code>include-path-deal2</code> </dt>
<dd> <p>
Same for the /deal.II library. Include files are in
libraries as <code>libg2c</code> for g77, or vendor libraries
</p>
</dd>
-
-
- <dt> <code>LIBS</code> </dt>
- <dd> <p>
- Additional libraries found by configure and needed to be
- linked with based on the flags given to
- <code>./configure</code>, e.g. BLAS or LAPACK
- </p>
- </dd>
</dl>
<h3>General</h3>
<ol>
+ <li> <p>
+ New: There is now some support to include and use routines from the
+ <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
+ target="_top">Harwell Subroutine Library</a>.
+ <br>
+ (WB 2001/01/30)
+ </p>
+
<li> <p>
New: The <code>./configure</code> script now checks for the
existence of a Fortran 77 compiler, and sets its path, some
<h3>lac</h3>
<ol>
+ <li> <p>
+ New: There is now some support to include and use routines from the
+ <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
+ target="_top">Harwell Subroutine Library</a>, and support
+ classes for the sparse direct solvers MA27 and MA47.
+ <br>
+ (WB 2001/01/30)
+ </p>
+
</ol>
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)
--- /dev/null
+//---------------------------- 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 <base/exceptions.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+
+
+
+/**
+ * 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<double> &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<double> &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<double> &matrix,
+ Vector<double> &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<const SparsityPattern> 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<unsigned int> row_numbers;
+ vector<unsigned int> column_numbers;
+
+ /**
+ * Array to hold the matrix
+ * elements, and later the
+ * elements of the factors.
+ */
+ vector<double> 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<unsigned int> IW;
+ vector<unsigned int> IKEEP;
+ vector<unsigned int> 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<double> &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<double> &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<double> &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<double> &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<double> &matrix,
+ Vector<double> &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<const SparseMatrix<double> > 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<unsigned int> row_numbers;
+ vector<unsigned int> column_numbers;
+
+ /**
+ * Array to hold the matrix
+ * elements, and later the
+ * elements of the factors.
+ */
+ vector<double> 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<unsigned int> IW;
+ vector<unsigned int> KEEP;
+ vector<unsigned int> IW1;
+
+ /**
+ * Fill the @p{A} array from the
+ * symmetric part of the given
+ * matrix.
+ */
+ void fill_A (const SparseMatrix<double> &matrix);
+};
+
+
+
+
+#endif
--- /dev/null
+//---------------------------- 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 <lac/sparse_direct.h>
+#include <hsl/hsl.h>
+#include <base/memory_consumption.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector.h>
+
+
+// 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<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ if (row <= *col)
+ ++n_nonzero_elements;
+
+
+ // fill the row numbers and column
+ // numbers arrays from the sparsity
+ // pattern. note that we have
+ // Fortran convention, i.e. indices
+ // need to be 1-base, as opposed to
+ // C's 0-based convention!
+ row_numbers.resize (n_nonzero_elements);
+ column_numbers.resize (n_nonzero_elements);
+
+ unsigned int global_index = 0;
+ for (unsigned int row=0; row<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ // note that the matrix must be
+ // symmetric, so only treat the
+ // upper right part
+ if (row <= *col)
+ {
+ Assert (global_index < n_nonzero_elements, ExcInternalError());
+
+ row_numbers[global_index] = row+1;
+ column_numbers[global_index] = *col+1;
+ ++global_index;
+ };
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
+
+ // initialize scratch arrays and
+ // variables
+ LIW = static_cast<unsigned int>((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<unsigned int>(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<double> &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<int>(NRLNEC * LA_factor);
+ A.resize (LA);
+ fill_A (matrix);
+
+ // if necessary extend IW
+ if (LIW < NIRNEC * LIW_factor_2)
+ {
+ LIW = static_cast<unsigned int>(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<unsigned int>(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<unsigned int>(LA * LA_increase_factor);
+ if (true)
+ {
+ vector<double> 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<double> &rhs_and_solution) const
+{
+ Assert (factorize_called == true, ExcFactorizeNotCalled());
+
+ const unsigned int n_rows = rhs_and_solution.size();
+ vector<double> 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<double> &matrix,
+ Vector<double> &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<double> &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<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ // note that the matrix must be
+ // symmetric, so only treat the
+ // upper right part
+ if (row <= *col)
+ {
+ Assert (global_index < n_nonzero_elements, ExcInternalError());
+
+ A[global_index] = matrix(row,*col);
+ ++global_index;
+ };
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
+};
+
+
+
+/* -------------------------- MA47 ---------------------------- */
+
+
+
+
+SparseDirectMA47::SparseDirectMA47 (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),
+ matrix (0)
+{
+ HSL::MA47::ma47id_ (CNTL, ICNTL);
+
+ // suppress error output if
+ // requested
+ if (suppress_output)
+ ICNTL[0] = 0;
+};
+
+
+
+void
+SparseDirectMA47::initialize (const SparseMatrix<double> &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<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ // skip zero elements, as
+ // required by the docs of MA47
+ if ((row <= *col) && (m(row,*col) != 0))
+ ++n_nonzero_elements;
+
+
+ // fill the row numbers and column
+ // numbers arrays from the sparsity
+ // pattern. note that we have
+ // Fortran convention, i.e. indices
+ // need to be 1-base, as opposed to
+ // C's 0-based convention!
+ row_numbers.resize (n_nonzero_elements);
+ column_numbers.resize (n_nonzero_elements);
+
+ unsigned int global_index = 0;
+ for (unsigned int row=0; row<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ // note that the matrix must be
+ // symmetric, so only treat the
+ // upper right part
+ if ((row <= *col) && (m(row,*col) != 0))
+ {
+ Assert (global_index < n_nonzero_elements, ExcInternalError());
+
+ row_numbers[global_index] = row+1;
+ column_numbers[global_index] = *col+1;
+ ++global_index;
+ };
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
+
+ // initialize scratch arrays and
+ // variables
+ LIW = static_cast<unsigned int>((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<unsigned int>(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<double> &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<int>(INFO[5] * LA_factor);
+ A.resize (LA);
+ fill_A (m);
+
+ // if necessary extend IW
+ if (LIW < INFO[6] * LIW_factor_2)
+ {
+ LIW = static_cast<unsigned int>(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<unsigned int>(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<unsigned int>(LA * LA_increase_factor);
+ if (true)
+ {
+ vector<double> 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<double> &rhs_and_solution)
+{
+ Assert (factorize_called == true, ExcFactorizeNotCalled());
+
+ const unsigned int n_rows = rhs_and_solution.size();
+ vector<double> 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<double> &matrix,
+ Vector<double> &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<double> &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<n_rows; ++row)
+ for (const unsigned int *col=&col_nums[rowstart_indices[row]];
+ col != &col_nums[rowstart_indices[row+1]];
+ ++col)
+ // note that the matrix must be
+ // symmetric, so only treat the
+ // upper right part
+ if ((row <= *col) && (matrix(row,*col) != 0))
+ {
+ Assert (global_index < n_nonzero_elements, ExcInternalError());
+
+ A[global_index] = matrix(row,*col);
+ ++global_index;
+ };
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
+};