]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Include support for the Harwell Subroutine Library (HSL) into the library, in particu...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Jan 2001 16:38:28 +0000 (16:38 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 30 Jan 2001 16:38:28 +0000 (16:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@3844 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
deal.II/Makefile
deal.II/common/Make.global_options.in
deal.II/configure
deal.II/configure.in
deal.II/contrib/Makefile.in
deal.II/contrib/configure.in
deal.II/contrib/hsl/lib/o/.cvsignore [new file with mode: 0644]
deal.II/doc/development/Makefile
deal.II/doc/development/makefiles.1.html
deal.II/doc/news/2001/c-3-1.html
deal.II/lac/Makefile
deal.II/lac/include/lac/sparse_direct.h [new file with mode: 0644]
deal.II/lac/source/sparse_direct.cc [new file with mode: 0644]

index 7f863cbed45a6c3ff1c54b67e8f86d99771995bb..2a116b96b2280815cf06664f2615e0a9a90ae0cd 100644 (file)
@@ -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
index 911c08df1aa199d2e6b19c08b8185af8ef0f28c2..1b5c0d16a36785802e8aeb9c8332fd14af0ea8b0 100644 (file)
@@ -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)
index 15a2f2a5c686ffc8c676cc45e61b0fcd94d0710e..31cb7e90a9e170f651b5c21e611e22b5cdb74e8e 100755 (executable)
@@ -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 <<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
@@ -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 <<EOF
-#line 2021 "configure"
+#line 2053 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 #ifdef __cplusplus
@@ -2031,7 +2063,7 @@ int main() {
 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
@@ -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 <<EOF
-#line 2196 "configure"
+#line 2228 "configure"
 #include "confdefs.h"
 
 #include <sys/resource.h>
@@ -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 <<EOF
-#line 2237 "configure"
+#line 2269 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2265 "configure"
+#line 2297 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2297 "configure"
+#line 2329 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2322 "configure"
+#line 2354 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2367 "configure"
+#line 2399 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2395 "configure"
+#line 2427 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2427 "configure"
+#line 2459 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<EOF
-#line 2452 "configure"
+#line 2484 "configure"
 #include "confdefs.h"
 
 #include <cmath>
@@ -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 <<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
@@ -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 <<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
@@ -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
index c1bfe5ede288583162e0185d6679d7ae083c2612..1e0867a60f27417c97cd7499effccee6605fce5d 100644 (file)
@@ -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!
index b70d5ab319fed8223df5d85d9e7fa0226641bb1c..9fa5eb3e2b0ec9ad60707e3dff6ea88a32342a78 100644 (file)
@@ -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
index f1f8759b322418606b4323d1fa810a480afe89ea..4e93042c349b8584e2c7a818b1e1d01c816cf859 100644 (file)
@@ -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/contrib/hsl/lib/o/.cvsignore b/deal.II/contrib/hsl/lib/o/.cvsignore
new file mode 100644 (file)
index 0000000..5761abc
--- /dev/null
@@ -0,0 +1 @@
+*.o
index 30f9c89ccb120b4ab42896c5b0dd8c2fb447e129..8002c8a7584a5c2c9d1f7025542a58522d50c97e 100644 (file)
@@ -19,9 +19,11 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
        @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>' >> $@
@@ -33,9 +35,11 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
        @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>' >> $@
@@ -43,7 +47,6 @@ makefiles.html: Makefile makefiles.1.html $D/common/Make.global_options makefile
        @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>' >> $@
index a8a5be3c26eb12b1d7df986b9f3dbac955b599bb..dd09c0351f18868a934da4ad880744938449057a 100644 (file)
@@ -85,7 +85,8 @@
            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>
 
 
index bee329ab9566c270c1f78ac87a8a1a1f7b9a8f39..4efeb7ce3f0fe6723ffce1b403464d1f0c0e5cc9 100644 (file)
@@ -26,6 +26,14 @@ documentation, etc</a>.
 <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
@@ -68,6 +76,15 @@ documentation, etc</a>.
 <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>
 
 
index 71d16407065355d73616a157fcdf4202d823201f..501f45a9d5559a7e276a15c487575d260794e5a6 100644 (file)
@@ -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 (file)
index 0000000..1330e7a
--- /dev/null
@@ -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 <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
diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc
new file mode 100644 (file)
index 0000000..b0879d5
--- /dev/null
@@ -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 <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());  
+};

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.