release: 0.1.1

26 files changed, 472 insertions, 205 deletions

diff --git a/flower/Makefile b/flower/Makefile
index 30b2a6f8f4..1089809e21 100644
--- a/flower/Makefile
+++ b/flower/Makefile
@@ -39,7 +39,7 @@ SUBDIRS = include test
 # list of distribution files:
 #
 SCRIPTS = 
-README_FILES = NEWS README TODO
+README_FILES = ONEWS NEWS README TODO
 EXTRA_DISTFILES= configure config.hh.in configure.in VERSION $(README_FILES) $(SCRIPTS) Flower-flags.make.in
 #
 
@@ -74,3 +74,6 @@ endif
 
 localuninstall:
 	rm -f $(libdir)/libflower.{so,a}
+
+
+$(outdir)/flower-version.o: $(outdir)/version.hh
diff --git a/flower/NEWS b/flower/NEWS
index a93f485302..57ad102788 100644
--- a/flower/NEWS
+++ b/flower/NEWS
@@ -1,4 +1,9 @@
-version 1.1:
+
+
+pl 25
+	- unordered_substitute and unordered_del
+	- Full_storage inline functions
+	- more Diagonal_storage fixes
 
 pl 24 
 	- Diagonal_storage for band matrices.
@@ -9,6 +14,7 @@ pl 24
 pl 23:
 	- virtual-methods : static_is_type_b
 	- Choleski checks off if not PARANOID
+
 pl 22:
 	- ACursor and PACursor to give array a List like interface.
 
@@ -18,6 +24,7 @@ pl 20.mb
 
 pl 20
 	- List::junk_links()
+
 pl 19
 	- Array::reverse()
 
@@ -29,6 +36,7 @@ pl 18
 
 pl 17
 	- naming: Pointer->Link, IPointer->Pointer
+
 pl 16
 	- Array::get()
 	- P< > doco.
@@ -39,6 +47,7 @@ pl 15
 
 pl 14
 	- interval methods
+
 pl 13
 	- better test-bed
 	- Heap PQueue implementation
@@ -115,36 +124,3 @@ pl 1
 0:
 		
 
-------------------
-1.0:
-
-pl 27-3
-	- debug memmove code
-	- StringData bugfix
-	- old String::String( char, int ) bugfix
-
-pl 27-1,2 (not released)
-	patches by JCN
-	- stringutils.cc included again 
-	- bin2hex_str bugfix
-	- String class handles null bytes
-	- StringUtils inlined/outlined by #define
-	- StringConversion (only hex for now)
-
-pl 27
-	- (temporarily?) removed findcursor* t*mat*
-
-pl 26
-	- docxx 3.0
-
-pl 25
-	- merge sstack and Array
-
-pl 24
-	- small fix in vector print
-
-pl 23
-	- win32 patches (JN)
-
-pl 22
-	- Array::add -> Array::push
diff --git a/flower/ONEWS b/flower/ONEWS
new file mode 100644
index 0000000000..9f4df5bb07
--- /dev/null
+++ b/flower/ONEWS
@@ -0,0 +1,34 @@
+history of flower lib 0.0
+
+version 1.1  release
+
+pl 27-3
+	- debug memmove code
+	- StringData bugfix
+	- old String::String( char, int ) bugfix
+
+pl 27-1,2 (not released)
+	patches by JCN
+	- stringutils.cc included again 
+	- bin2hex_str bugfix
+	- String class handles null bytes
+	- StringUtils inlined/outlined by #define
+	- StringConversion (only hex for now)
+
+pl 27
+	- (temporarily?) removed findcursor* t*mat*
+
+pl 26
+	- docxx 3.0
+
+pl 25
+	- merge sstack and Array
+
+pl 24
+	- small fix in vector print
+
+pl 23
+	- win32 patches (JN)
+
+pl 22
+	- Array::add -> Array::push
diff --git a/flower/TODO b/flower/TODO
index 169a0525e6..b71ad836c0 100644
--- a/flower/TODO
+++ b/flower/TODO
@@ -1,5 +1,8 @@
 	* write a String_hash template
 
+
+	* Array::slice() upper too
+
 	* write a Pointer_hash template
 
 	* fix/junk ambiguous String constructor overloads, e.g.:
diff --git a/flower/VERSION b/flower/VERSION
index fafe5106a4..d05f7ff918 100644
--- a/flower/VERSION
+++ b/flower/VERSION
@@ -1,6 +1,6 @@
 MAJOR_VERSION = 1
 MINOR_VERSION = 1
-PATCH_LEVEL = 24
+PATCH_LEVEL = 25
 # use to send patches, always empty for released version:
 MY_PATCH_LEVEL = # include separator: "-1" or ".a"
 #
diff --git a/flower/choleski.cc b/flower/choleski.cc
index d619b97460..086759193f 100644
--- a/flower/choleski.cc
+++ b/flower/choleski.cc
@@ -12,13 +12,15 @@ const Real EPS = 1e-7;		// so sue me. Hard coded
 // for testing new Matrix_storage.
 //#define PARANOID
 
-Vector
-Choleski_decomposition::solve(Vector rhs)const
+void
+Choleski_decomposition::full_matrix_solve(Vector &out, Vector const &rhs)const
 {
     int n= rhs.dim();
     assert(n == L.dim());
-    Vector y(n);
-
+    Vector y;
+    y.set_dim( n);
+    out.set_dim(n);
+    
     // forward substitution
     for (int i=0; i < n; i++) {
 	Real sum(0.0);
@@ -26,18 +28,67 @@ Choleski_decomposition::solve(Vector rhs)const
 	    sum += y(j) * L(i,j);
 	y(i) = (rhs(i) - sum)/L(i,i);
     }
+    
     for (int i=0; i < n; i++)
 	y(i) /= D(i);
 
     // backward subst
-    Vector &x(rhs);		// using input as return val.
+    Vector &x(out);		// using input as return val.
     for (int i=n-1; i >= 0; i--) {
 	Real sum(0.0);
 	for (int j=i+1; j < n; j++)
 	    sum += L(j,i)*x(j);
 	x(i) = (y(i) - sum)/L(i,i);
     }
-    return x;
+}
+
+void
+Choleski_decomposition::band_matrix_solve(Vector &out, Vector const &rhs)const
+{
+    int n= rhs.dim();
+    int b = L.band_i();
+    assert(n == L.dim());
+
+      out.set_dim(n);
+  
+    Vector y;
+    y.set_dim(n);
+    
+    // forward substitution
+    for (int i=0; i < n; i++) {
+	Real sum(0.0);
+	for (int j= 0 >? i - b; j < i; j++)
+	    sum += y(j) * L(i,j);
+	y(i) = (rhs(i) - sum)/L(i,i);
+    }
+    for (int i=0; i < n; i++)
+	y(i) /= D(i);
+
+    // backward subst
+    Vector &x(out);		// using input as return val.
+    for (int i=n-1; i >= 0; i--) {
+	Real sum(0.0);
+	for (int j=i+1; j <= i + b&&j < n ; j++)
+	    sum += L(j,i)*x(j);
+	x(i) = (y(i) - sum)/L(i,i);
+    }
+}
+
+void
+Choleski_decomposition::solve(Vector &x, Vector const &rhs)const
+{
+    if (L.band_b()) {
+	band_matrix_solve(x,rhs);
+    } else
+	full_matrix_solve(x,rhs);
+}
+
+Vector
+Choleski_decomposition::solve(Vector rhs)const
+{
+    Vector r;
+    solve(r, rhs);
+    return r;
 }
 
 void
@@ -126,9 +177,10 @@ Choleski_decomposition::inverse() const
     int n=L.dim();
     Matrix invm(n);
     Vector e_i(n);
+    Vector inv(n);
     for (int i = 0; i < n; i++) {
 	e_i.set_unit(i);
-	Vector inv(solve(e_i));
+	solve(inv, e_i);
 	for (int j = 0 ; j<n; j++)
 	    invm(i,j) = inv(j);
     }
diff --git a/flower/configure b/flower/configure
index 53ff4301cf..920d8c0c08 100755
--- a/flower/configure
+++ b/flower/configure
@@ -15,6 +15,8 @@ ac_help="$ac_help
   enable-shared       shared flower library"
 ac_help="$ac_help
   disable-optimise    optimisations off"
+ac_help="$ac_help
+  out-dir			set the directory for machine generated files. Default out or out-HOST"
 
 # Initialize some variables set by options.
 # The variables have the same names as the options, with
@@ -534,16 +536,18 @@ optimise_b=yes
 shared_b=no
 LIB_SUFFIX=.a
 
-# if given here, these vars are initted at the checking point.
 if test x$host = xNONE; then
-    flowerbuildprefix=.
+    OUTDIR_NAME=${OUTDIR_NAME-"out"}
 else
-    flowerbuildprefix="../$host-build-dir/Flower"
-    mkdir $flowerbuildprefix;
-    for a in `find -type d -and -not -name '*-build-dir'`; do
-	mkdir $flowerbuildprefix/$a;
-    done
+    OUTDIR_NAME=${OUTDIR_NAME-"out-$host"}
 fi
+
+for a in `find -type d -and -not -name 'out'`; do
+	if test ! -d  $a/$OUTDIR_NAME; then
+		mkdir $a/$OUTDIR_NAME;
+	fi
+done
+
     
 # Check whether --enable-shared or --disable-shared was given.
 if test "${enable_shared+set}" = set; then
@@ -559,6 +563,14 @@ if test "${enable_optimise+set}" = set; then
 fi
 
 
+# Check whether --enable-out-dir or --disable-out-dir was given.
+if test "${enable_out_dir+set}" = set; then
+  enableval="$enable_out_dir"
+  OUTDIR_NAME=$enableval
+    
+fi
+
+
 if test $shared_b = yes; then
 	MODULE_CXXFLAGS="$MODULE_CXXFLAGS -fPIC"
 	MODULE_LDFLAGS="-shared -Wl,-soname,libflower.so "
@@ -578,7 +590,7 @@ do
 # Extract the first word of "$ac_prog", so it can be a program name with args.
 set dummy $ac_prog; ac_word=$2
 echo $ac_n "checking for $ac_word""... $ac_c" 1>&6
-echo "configure:582: checking for $ac_word" >&5
+echo "configure:594: checking for $ac_word" >&5
 if eval "test \"`echo '$''{'ac_cv_prog_CXX'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
@@ -609,7 +621,7 @@ test -n "$CXX" || CXX="gcc"
 
 
 echo $ac_n "checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) works""... $ac_c" 1>&6
-echo "configure:613: checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) works" >&5
+echo "configure:625: checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) works" >&5
 
 ac_ext=C
 # CXXFLAGS is not in ac_cpp because -g, -O, etc. are not valid cpp options.
@@ -619,11 +631,11 @@ ac_link='${CXX-g++} -o conftest $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $L
 cross_compiling=$ac_cv_prog_cxx_cross
 
 cat > conftest.$ac_ext <<EOF
-#line 623 "configure"
+#line 635 "configure"
 #include "confdefs.h"
 main(){return(0);}
 EOF
-if { (eval echo configure:627: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
+if { (eval echo configure:639: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
   ac_cv_prog_cxx_works=yes
   # If we can't run a trivial program, we are probably using a cross compiler.
   if (./conftest; exit) 2>/dev/null; then
@@ -643,12 +655,12 @@ if test $ac_cv_prog_cxx_works = no; then
   { echo "configure: error: installation or configuration problem: C++ compiler cannot create executables." 1>&2; exit 1; }
 fi
 echo $ac_n "checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) is a cross-compiler""... $ac_c" 1>&6
-echo "configure:647: checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) is a cross-compiler" >&5
+echo "configure:659: checking whether the C++ compiler ($CXX $CXXFLAGS $LDFLAGS) is a cross-compiler" >&5
 echo "$ac_t""$ac_cv_prog_cxx_cross" 1>&6
 cross_compiling=$ac_cv_prog_cxx_cross
 
 echo $ac_n "checking whether we are using GNU C++""... $ac_c" 1>&6
-echo "configure:652: checking whether we are using GNU C++" >&5
+echo "configure:664: checking whether we are using GNU C++" >&5
 if eval "test \"`echo '$''{'ac_cv_prog_gxx'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
@@ -657,7 +669,7 @@ else
   yes;
 #endif
 EOF
-if { ac_try='${CXX-g++} -E conftest.C'; { (eval echo configure:661: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; }; } | egrep yes >/dev/null 2>&1; then
+if { ac_try='${CXX-g++} -E conftest.C'; { (eval echo configure:673: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; }; } | egrep yes >/dev/null 2>&1; then
   ac_cv_prog_gxx=yes
 else
   ac_cv_prog_gxx=no
@@ -672,7 +684,7 @@ if test $ac_cv_prog_gxx = yes; then
   ac_save_CXXFLAGS="$CXXFLAGS"
   CXXFLAGS=
   echo $ac_n "checking whether ${CXX-g++} accepts -g""... $ac_c" 1>&6
-echo "configure:676: checking whether ${CXX-g++} accepts -g" >&5
+echo "configure:688: checking whether ${CXX-g++} accepts -g" >&5
 if eval "test \"`echo '$''{'ac_cv_prog_cxx_g'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
@@ -700,7 +712,7 @@ else
 fi
 
 echo $ac_n "checking for 8-bit clean memcmp""... $ac_c" 1>&6
-echo "configure:704: checking for 8-bit clean memcmp" >&5
+echo "configure:716: checking for 8-bit clean memcmp" >&5
 if eval "test \"`echo '$''{'ac_cv_func_memcmp_clean'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
@@ -708,7 +720,7 @@ else
   ac_cv_func_memcmp_clean=no
 else
   cat > conftest.$ac_ext <<EOF
-#line 712 "configure"
+#line 724 "configure"
 #include "confdefs.h"
 #ifdef __cplusplus
 extern "C" void exit(int);
@@ -721,7 +733,7 @@ main()
 }
 
 EOF
-if { (eval echo configure:725: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest && (./conftest; exit) 2>/dev/null
+if { (eval echo configure:737: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest && (./conftest; exit) 2>/dev/null
 then
   ac_cv_func_memcmp_clean=yes
 else
@@ -739,12 +751,12 @@ echo "$ac_t""$ac_cv_func_memcmp_clean" 1>&6
 test $ac_cv_func_memcmp_clean = no && LIBOBJS="$LIBOBJS memcmp.o"
 
 echo $ac_n "checking for vprintf""... $ac_c" 1>&6
-echo "configure:743: checking for vprintf" >&5
+echo "configure:755: checking for vprintf" >&5
 if eval "test \"`echo '$''{'ac_cv_func_vprintf'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 748 "configure"
+#line 760 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char vprintf(); below.  */
@@ -770,7 +782,7 @@ vprintf();
 
 ; return 0; }
 EOF
-if { (eval echo configure:774: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
+if { (eval echo configure:786: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
   rm -rf conftest*
   eval "ac_cv_func_vprintf=yes"
 else
@@ -794,12 +806,12 @@ fi
 
 if test "$ac_cv_func_vprintf" != yes; then
 echo $ac_n "checking for _doprnt""... $ac_c" 1>&6
-echo "configure:798: checking for _doprnt" >&5
+echo "configure:810: checking for _doprnt" >&5
 if eval "test \"`echo '$''{'ac_cv_func__doprnt'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 803 "configure"
+#line 815 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char _doprnt(); below.  */
@@ -825,7 +837,7 @@ _doprnt();
 
 ; return 0; }
 EOF
-if { (eval echo configure:829: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
+if { (eval echo configure:841: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
   rm -rf conftest*
   eval "ac_cv_func__doprnt=yes"
 else
@@ -852,12 +864,12 @@ fi
 for ac_func in memmem snprintf 
 do
 echo $ac_n "checking for $ac_func""... $ac_c" 1>&6
-echo "configure:856: checking for $ac_func" >&5
+echo "configure:868: checking for $ac_func" >&5
 if eval "test \"`echo '$''{'ac_cv_func_$ac_func'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 861 "configure"
+#line 873 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char $ac_func(); below.  */
@@ -883,7 +895,7 @@ $ac_func();
 
 ; return 0; }
 EOF
-if { (eval echo configure:887: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
+if { (eval echo configure:899: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest; then
   rm -rf conftest*
   eval "ac_cv_func_$ac_func=yes"
 else
@@ -908,8 +920,10 @@ fi
 done
 
 
+if test ! -d $OUTDIR_NAME ; then
+    mkdir $OUTDIR_NAME
+fi
 
-CXX="$ac_cv_prog_CXX" ../bin/make-version > $flowerbuildprefix/out/version.hh
 
 trap '' 1 2 15
 cat > confcache <<\EOF
@@ -1011,7 +1025,7 @@ done
 
 ac_given_srcdir=$srcdir
 
-trap 'rm -fr `echo "$flowerbuildprefix/out/Flower-flags.make:Flower-flags.make.in $flowerbuildprefix/out/config.hh:config.hh.in" | sed "s/:[^ ]*//g"` conftest*; exit 1' 1 2 15
+trap 'rm -fr `echo "$OUTDIR_NAME/Flower-flags.make:Flower-flags.make.in $OUTDIR_NAME/config.hh:config.hh.in" | sed "s/:[^ ]*//g"` conftest*; exit 1' 1 2 15
 EOF
 cat >> $CONFIG_STATUS <<EOF
 
@@ -1087,7 +1101,7 @@ EOF
 
 cat >> $CONFIG_STATUS <<EOF
 
-CONFIG_FILES=\${CONFIG_FILES-"$flowerbuildprefix/out/Flower-flags.make:Flower-flags.make.in"}
+CONFIG_FILES=\${CONFIG_FILES-"$OUTDIR_NAME/Flower-flags.make:Flower-flags.make.in"}
 EOF
 cat >> $CONFIG_STATUS <<\EOF
 for ac_file in .. $CONFIG_FILES; do if test "x$ac_file" != x..; then
@@ -1163,7 +1177,7 @@ ac_eD='%g'
 if test "${CONFIG_HEADERS+set}" != set; then
 EOF
 cat >> $CONFIG_STATUS <<EOF
-  CONFIG_HEADERS="$flowerbuildprefix/out/config.hh:config.hh.in"
+  CONFIG_HEADERS="$OUTDIR_NAME/config.hh:config.hh.in"
 EOF
 cat >> $CONFIG_STATUS <<\EOF
 fi
diff --git a/flower/configure.in b/flower/configure.in
index 619a5fe71e..dd20827eae 100644
--- a/flower/configure.in
+++ b/flower/configure.in
@@ -7,16 +7,18 @@ optimise_b=yes
 shared_b=no
 LIB_SUFFIX=.a
 
-# if given here, these vars are initted at the checking point.
 if test x$host = xNONE; then
-    flowerbuildprefix=.
+    OUTDIR_NAME=${OUTDIR_NAME-"out"}
 else
-    flowerbuildprefix="../$host-build-dir/Flower"
-    mkdir $flowerbuildprefix;
-    for a in `find -type d -and -not -name '*-build-dir'`; do
-	mkdir $flowerbuildprefix/$a;
-    done
+    OUTDIR_NAME=${OUTDIR_NAME-"out-$host"}
 fi
+
+for a in `find -type d -and -not -name 'out'`; do
+	if test ! -d  $a/$OUTDIR_NAME; then
+		mkdir $a/$OUTDIR_NAME;
+	fi
+done
+
     
 AC_ARG_ENABLE(shared,
     [  enable-shared       shared flower library],
@@ -26,6 +28,11 @@ AC_ARG_ENABLE(optimise,
     [  disable-optimise    optimisations off],
     [optimise_b=$enableval])
 
+AC_ARG_ENABLE(out-dir,
+    [  out-dir			set the directory for machine generated files. Default out or out-HOST],
+    [OUTDIR_NAME=$enableval]
+    [])
+
 if test $shared_b = yes; then
 	MODULE_CXXFLAGS="$MODULE_CXXFLAGS -fPIC"
 	MODULE_LDFLAGS="-shared -Wl,-soname,libflower.so "
@@ -44,9 +51,11 @@ AC_PROG_CXX
 AC_FUNC_MEMCMP
 AC_FUNC_VPRINTF
 AC_CHECK_FUNCS(memmem snprintf )
-AC_CONFIG_HEADER($flowerbuildprefix/out/config.hh:config.hh.in)
 
-CXX="$ac_cv_prog_CXX" ../bin/make-version > $flowerbuildprefix/out/version.hh
+if test ! -d $OUTDIR_NAME ; then
+    mkdir $OUTDIR_NAME
+fi
 
-AC_OUTPUT($flowerbuildprefix/out/Flower-flags.make:Flower-flags.make.in)
+AC_CONFIG_HEADER($OUTDIR_NAME/config.hh:config.hh.in)
+AC_OUTPUT($OUTDIR_NAME/Flower-flags.make:Flower-flags.make.in)
 
diff --git a/flower/diagonal-storage.cc b/flower/diagonal-storage.cc
index 3cf2f93d20..fc812423e5 100644
--- a/flower/diagonal-storage.cc
+++ b/flower/diagonal-storage.cc
@@ -9,6 +9,15 @@
 
 #include "diagonal-storage.hh"
 
+
+#ifdef INLINE
+#undef INLINE
+#endif
+
+#define INLINE inline
+
+#include "full-storage.icc"
+
 int
 Diagonal_storage::dim()const
 {
@@ -40,6 +49,7 @@ Diagonal_storage::band_size_i()const
 void
 Diagonal_storage::set_band_size(int s)
 {
+    assert( s>=0);
     Full_storage f(dim(), 2*s+1);
     for (int i=0; i < dim(); i++) {
 	int k=-s;
@@ -110,7 +120,7 @@ void
 Diagonal_storage::resize_dim(int d)
 {
     Full_storage f(d, 2*band_size_i()+1);
-    for (int i=0; i < d&& i < dim(); i++) {
+    for (int i=0; i < d && i < dim(); i++) {
 	for ( int k=0;  k < 2*band_size_i(); k++)
 	    f.elem(i,k) = elem(i,k);
     }
@@ -134,9 +144,7 @@ Diagonal_storage::mult_next(int &i, int &j)const
 	j = i- band_size_i();
     if ( j > i + band_size_i() || j >= dim() ) {
 	i++;
-	j = i - band_size_i(); 
-	if  (j < 0)
-	    j=0;
+	j = 0 >? i - band_size_i(); 
     }
 }
 
@@ -155,9 +163,7 @@ Diagonal_storage::trans_next(int &i, int& j)const
     
     if ( i >= dim() || i > j + band_size_i() ) {
 	j++;
-	i = j - band_size_i(); 
-	if  (i < 0)
-	    i=0;
+	i = 0 >? j - band_size_i(); 
     }
 }
 
@@ -178,8 +184,9 @@ Diagonal_storage::elem(int i, int j)
     /*
       if this fails, the previous call fucked up
       */
-    assert(nul_entry);
-    if (abs ( i-j ) > band_size_i())
+    assert(!nul_entry);
+
+    if (abs ( i-j ) > band_size_i())  
 	return nul_entry;
     else
 	return band_.elem(i, j - i + band_size_i());
@@ -196,8 +203,8 @@ Diagonal_storage::try_right_multiply(Matrix_storage*dest,
     if ( right->name() != Diagonal_storage::static_name() ) 
 	return false;
     
-    const Diagonal_storage*  diag = (Diagonal_storage const*)right;
-    int band2 = diag->band_size_i();
+    const Diagonal_storage*  right_diag = (Diagonal_storage const*)right;
+    int band2 = right_diag->band_size_i();
     int n = dim();
     /*
       should check if dest is a Diagonal_storage of sufficient size too.
@@ -209,7 +216,7 @@ Diagonal_storage::try_right_multiply(Matrix_storage*dest,
 	    int relk =  startk + band_size_i() -i;
 	    Real sum =0.0;
 	    for ( int k = startk; k <= stopk; k++)
-		sum += band_.elem(i, relk) * diag->elem(relk, j);
+		sum += band_.elem(i, relk++) * right_diag->elem(k, j);
 	    dest->elem(i, j) = sum;
 	    
 	}
@@ -234,3 +241,5 @@ Diagonal_storage::OK() const
 {
     band_.OK();
 }
+
+IMPLEMENT_VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);
diff --git a/flower/directed-graph.cc b/flower/directed-graph.cc
index c47285c71d..a08f456381 100644
--- a/flower/directed-graph.cc
+++ b/flower/directed-graph.cc
@@ -68,7 +68,7 @@ Directed_graph_node::remove_edge_out_idx(int i)
 
     int j = d_l->edge_in_l_arr_.find_i(this);
     assert(j>=0);
-    d_l->edge_in_l_arr_.del(j);
+    d_l->edge_in_l_arr_.unordered_del(j);
     PARANOID_OK();
 }
 
diff --git a/flower/full-storage.cc b/flower/full-storage.cc
index 223796ba35..fd0d9f550b 100644
--- a/flower/full-storage.cc
+++ b/flower/full-storage.cc
@@ -8,6 +8,7 @@
 
 #include "full-storage.hh"
 
+
 void
 Full_storage::operator=(Full_storage const &fs)
 {
@@ -19,6 +20,7 @@ Full_storage::operator=(Full_storage const &fs)
 	    els_p_p_[i][j]= fs.els_p_p_[i][j];
 }
 
+
 void
 Full_storage::OK() const
 {
@@ -30,63 +32,18 @@ Full_storage::OK() const
 #endif
 }
 
-void
-Full_storage::resize_cols(int newh)
-{
-    if (newh <= max_height_i_) {
-	height_i_=newh;
-	return;
-    }
-   
-    Real ** newa=new Real*[newh];
-    int j=0;
-    for (; j < height_i_; j++)
-	newa[j] = els_p_p_[j];
-    for (; j < newh; j++)
-	newa[j] = new Real[max_width_i_];
-    delete[] els_p_p_;
-    els_p_p_=newa;
 
-    height_i_ = max_height_i_ = newh;
-}
 
 
-Full_storage::Full_storage(Matrix_storage*m)
+Full_storage::~Full_storage() 
 {
-    set_size(m->rows(), m->cols());
-    if ( !m->is_type_b ( Full_storage::static_name()))
-	for (int i=0; i<height_i_; i++)
-	    for (int j=0; j<width_i_; j++)
-		els_p_p_[i][j]=0.0;
-    for (int i,j=0; m->mult_ok(i,j); m->mult_next(i,j))
-	els_p_p_[i][j] = m->elem(i,j);
-}
-
-void
-Full_storage::resize_rows(int neww)
-{
-    if (neww <= max_width_i_) {
-	width_i_=neww;
-	return;
-    }
-    for (int i=0; i < max_height_i_ ; i++) {
-	Real* newa = new Real[neww];
-	for (int k=0; k < width_i_; k++)
-	    newa[k] = els_p_p_[i][k];
-
-	delete[] els_p_p_[i];
-	els_p_p_[i] = newa;
-    }
-    width_i_ = max_width_i_ = neww;	
-}
-
-Full_storage::~Full_storage() {
     for (int i=0; i < max_height_i_; i++)
 	delete [] els_p_p_[i];
     delete[] els_p_p_;
 }
 
 void
+
 Full_storage::resize(int rows, int cols)
 {
     OK();
@@ -95,12 +52,14 @@ Full_storage::resize(int rows, int cols)
 }
 
 
+
 bool
 Full_storage::mult_ok(int i, int ) const
 {
     return i < height_i_;
 }
 
+
 bool
 Full_storage::trans_ok(int , int j) const
 {
@@ -108,6 +67,7 @@ Full_storage::trans_ok(int , int j) const
 } 
 
 
+
 void
 Full_storage::trans_next(int &i, int &j) const
 {
@@ -119,6 +79,7 @@ Full_storage::trans_next(int &i, int &j) const
     }
 }
 
+
 void
 Full_storage::mult_next(int &i, int &j) const
 {
@@ -130,6 +91,7 @@ Full_storage::mult_next(int &i, int &j) const
     }
 }
 
+
 void
 Full_storage::delete_column(int k)
 {
@@ -139,6 +101,8 @@ Full_storage::delete_column(int k)
 	    els_p_p_[i][j-1]=els_p_p_[i][j];
     width_i_--;
 }
+
+
 void
 Full_storage::delete_row(int k)
 {
@@ -150,6 +114,7 @@ Full_storage::delete_row(int k)
 }
 
 
+
 void
 Full_storage::insert_row(int k)
 {
@@ -161,19 +126,6 @@ Full_storage::insert_row(int k)
 
 }
 
-int
-Full_storage::dim()const
-{
-    assert (rows()==cols());
-    return rows();
-}
-
-Full_storage::Full_storage(Full_storage const&s)
-{
-    init();
-    (*this) = s;
-}
-
 bool
 Full_storage::try_right_multiply(Matrix_storage * dest, Matrix_storage const * right)const
 {
@@ -198,3 +150,61 @@ Full_storage::try_right_multiply(Matrix_storage * dest, Matrix_storage const * r
     
 }
 IMPLEMENT_IS_TYPE_B1(Full_storage,Matrix_storage);
+void
+Full_storage::resize_cols(int newh)
+{
+    if (newh <= max_height_i_) {
+	height_i_=newh;
+	return;
+    }
+   
+    Real ** newa=new Real*[newh];
+    int j=0;
+    for (; j < height_i_; j++)
+	newa[j] = els_p_p_[j];
+    for (; j < newh; j++)
+	newa[j] = new Real[max_width_i_];
+    delete[] els_p_p_;
+    els_p_p_=newa;
+
+    height_i_ = max_height_i_ = newh;
+}
+
+
+
+Full_storage::Full_storage(Matrix_storage*m)
+{
+    set_size(m->rows(), m->cols());
+    if ( !m->is_type_b ( Full_storage::static_name()))
+	for (int i=0; i<height_i_; i++)
+	    for (int j=0; j<width_i_; j++)
+		els_p_p_[i][j]=0.0;
+    for (int i,j=0; m->mult_ok(i,j); m->mult_next(i,j))
+	els_p_p_[i][j] = m->elem(i,j);
+}
+
+
+void
+Full_storage::resize_rows(int neww)
+{
+    if (neww <= max_width_i_) {
+	width_i_=neww;
+	return;
+    }
+    for (int i=0; i < max_height_i_ ; i++) {
+	Real* newa = new Real[neww];
+	for (int k=0; k < width_i_; k++)
+	    newa[k] = els_p_p_[i][k];
+
+	delete[] els_p_p_[i];
+	els_p_p_[i] = newa;
+    }
+    width_i_ = max_width_i_ = neww;	
+}
+
+#ifdef INLINE
+#undef INLINE
+#endif
+#define INLINE
+
+#include "full-storage.icc"
diff --git a/flower/include/Makefile b/flower/include/Makefile
index 1ad45cb463..b7716f2676 100644
--- a/flower/include/Makefile
+++ b/flower/include/Makefile
@@ -14,5 +14,6 @@ include ./$(depth)/flower/VERSION
 #
 MODULE_NAME = flower
 
+
 #
 
diff --git a/flower/include/choleski.hh b/flower/include/choleski.hh
index a3a02a11e0..abd19ee676 100644
--- a/flower/include/choleski.hh
+++ b/flower/include/choleski.hh
@@ -30,7 +30,7 @@ struct Choleski_decomposition {
     solve Px = rhs
     */
     Vector solve(Vector rhs) const;
-
+    void solve (Vector &dest, Vector const &rhs)const;
     Vector operator * (Vector rhs) const { return solve (rhs); }
   /**
     return the inverse of the matrix P.
@@ -41,6 +41,8 @@ struct Choleski_decomposition {
     */
     Matrix original() const;
 private:
+    void full_matrix_solve(Vector &,Vector const&)const;
+    void band_matrix_solve(Vector &, Vector const&)const;
     void full_matrix_decompose(Matrix const & P);
     void band_matrix_decompose(Matrix const &P);
 	 
diff --git a/flower/include/diagonal-storage.hh b/flower/include/diagonal-storage.hh
index 517c5cd814..f4515edfb1 100644
--- a/flower/include/diagonal-storage.hh
+++ b/flower/include/diagonal-storage.hh
@@ -12,9 +12,9 @@
 #include "full-storage.hh"
 
 /**
-  Store a single-band matrix;
+  Store a  matrix with a single-band.
 
-  INVARIANT
+  @invariant
   
   Diagonal_storage(i,j) == band_(i, j-i + band_size_i())
   
@@ -58,7 +58,7 @@ public:
     virtual void mult_next(int &i, int &j) const ;
     virtual bool trans_ok(int i, int j) const;
     virtual void trans_next(int &i, int &j) const;
-    VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);
+    DECLARE_VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);
     DECLARE_MY_RUNTIME_TYPEINFO;
     virtual bool try_right_multiply(Matrix_storage * dest, Matrix_storage const *)const;
 };
diff --git a/flower/include/full-storage.hh b/flower/include/full-storage.hh
index 78bafab43d..1dfb96e8fa 100644
--- a/flower/include/full-storage.hh
+++ b/flower/include/full-storage.hh
@@ -26,57 +26,28 @@ class Full_storage : public Matrix_storage
     Real** els_p_p_;
 
     void
-    init() {
-	els_p_p_=0;
-	height_i_=width_i_=max_height_i_=max_width_i_=0;
-
-    }
-    
-    bool valid(int i, int j) const {
-	return (i>=0 && i < height_i_)
-	    && (j < width_i_ && j >=0);
-    }
-    
+    init() ;
+     
+    bool valid(int i, int j) const ; 
 
     void resize_rows(int);
     void resize_cols(int);
 
 public:
-    virtual int rows() const {
-	return height_i_;
-    }
-    virtual int cols() const {
-	return width_i_;
-    }
+    virtual int rows() const;
+    virtual int cols() const ;
     
     
     virtual void resize(int i, int j);
-    virtual void resize(int i) {
-	resize(i,i);
-    }
-    
-    virtual Real& elem(int i,int j) {
-	assert(valid(i,j));
-	return els_p_p_[i][j];
-    }
-    virtual Real elem(int i, int j) const {
-	assert(valid(i,j));
-	return els_p_p_[i][j];
-    }
+    virtual void resize(int i);
+    virtual Real& elem(int i,int j);
+    virtual Real elem(int i, int j)const ;
     int dim()const;
     Full_storage(Matrix_storage*);
-    Full_storage() {
-	init();
-    }
-    Full_storage(int i, int j) {
-	init();
-	set_size(i,j);
-    }
+    Full_storage();
+    Full_storage(int i, int j);
     Full_storage(Full_storage const&);
-    Full_storage(int i) {
-	init();
-	set_size(i);
-    }
+    Full_storage(int i);
     void OK() const;
     void operator=(Full_storage const &);
     
@@ -89,7 +60,7 @@ public:
     virtual void mult_next(int &i, int &j) const ;
     virtual bool trans_ok(int i, int j) const;
     virtual void trans_next(int &i, int &j) const;
-    VIRTUAL_COPY_CONS(Full_storage,Matrix_storage);
+    DECLARE_VIRTUAL_COPY_CONS(Full_storage,Matrix_storage);
     DECLARE_MY_RUNTIME_TYPEINFO;
     virtual bool try_right_multiply(Matrix_storage * dest, Matrix_storage const * )const;
 };
diff --git a/flower/include/full-storage.icc b/flower/include/full-storage.icc
new file mode 100644
index 0000000000..06bcc9d70c
--- /dev/null
+++ b/flower/include/full-storage.icc
@@ -0,0 +1,94 @@
+/*
+  full-storage.icc -- implement Full_storage inline functions
+
+  source file of the Flower Library
+
+  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+
+#ifndef FULL_STORAGE_ICC
+#define FULL_STORAGE_ICC
+
+
+INLINE void
+Full_storage::init()
+{
+    els_p_p_=0;
+    height_i_=width_i_=max_height_i_=max_width_i_=0;
+}
+INLINE bool
+Full_storage::valid(int i, int j)const
+{
+    return (i>=0 && i < height_i_)
+	&& (j < width_i_ && j >=0);
+}
+
+
+INLINE
+Full_storage::Full_storage(Full_storage const&s)
+{
+    init();
+    (*this) = s;
+}
+
+INLINE Real& 
+Full_storage::elem(int i,int j)
+{
+    assert(valid(i,j));
+    return els_p_p_[i][j];
+}
+
+INLINE Real
+Full_storage::elem(int i, int j) const {
+    assert(valid(i,j));
+    return els_p_p_[i][j];
+}
+
+INLINE
+Full_storage::Full_storage() {
+    init();
+}
+
+
+INLINE int
+Full_storage::rows() const
+{
+    return height_i_;
+}
+INLINE int
+Full_storage::cols() const
+{
+    return width_i_;
+}
+INLINE int
+Full_storage::dim()const
+{
+    assert (rows()==cols());
+    return rows();
+}
+
+INLINE void
+Full_storage::resize(int i)
+{
+    resize(i,i);
+}
+
+INLINE   
+Full_storage::Full_storage(int i,int j)
+{
+    init();
+    set_size(i,j);
+}
+
+INLINE 
+Full_storage::Full_storage(int i)
+{
+    init();
+    set_size(i);
+}
+ 
+INLINE
+IMPLEMENT_VIRTUAL_COPY_CONS(Full_storage,Matrix_storage);
+
+#endif // FULL_STORAGE_ICC
diff --git a/flower/include/parray.hh b/flower/include/parray.hh
index 1e7404ee55..61c69495c2 100644
--- a/flower/include/parray.hh
+++ b/flower/include/parray.hh
@@ -39,6 +39,17 @@ public:
 	    else
 		del(i);
     }
+    void unordered_substitute(T* old, T * new_l)
+    {
+	int i;
+	while ((i = find_i(old)) >=0) 
+	    if (new_l)
+		elem(i) =new_l;
+	    else {
+		unordered_del( i );
+	    }
+    
+    }
     void default_sort() {
 	sort(default_compare);
     }
diff --git a/flower/include/rational.hh b/flower/include/rational.hh
index 30fa2cb6d6..d3b68f726f 100644
--- a/flower/include/rational.hh
+++ b/flower/include/rational.hh
@@ -1 +1,17 @@
+/*
+  rational.hh -- declare 
+
+  source file of the Flower Library
+
+  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+
+#ifndef RATIONAL_HH
+#define RATIONAL_HH
 #include <Rational.h>
+
+/// print a Rational. To be called from the  debugger
+void print_rat(Rational const&);
+
+#endif // RATIONAL_HH
diff --git a/flower/include/real.hh b/flower/include/real.hh
index 0533af2111..457b82f967 100644
--- a/flower/include/real.hh
+++ b/flower/include/real.hh
@@ -1,3 +1,12 @@
+/*
+  real.hh -- declare Real
+
+  source file of the Flower Library
+
+  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+
 #ifndef REAL_HH
 #define REAL_HH
 
@@ -5,12 +14,15 @@
 #include <builtin.h>
 #include <minmax.h>
 #include <math.h>
+#include <limits.h>
 
 typedef double Real;
+const Real infinity_f = HUGE_VAL;
 
 inline Real
 distance(Real x,Real y)
 {
     return abs(x-y);
 }
+
 #endif
diff --git a/flower/include/varray.hh b/flower/include/varray.hh
index 5fe5fa3943..6f0e18aedc 100644
--- a/flower/include/varray.hh
+++ b/flower/include/varray.hh
@@ -164,6 +164,11 @@ public:
 	del (i);
 	return t;
     }
+    void unordered_del(int i)
+    {
+	elem(i) = top();
+	set_size(size() -1);
+    }
     void del(int i) {
 	assert(i >=0&& i < size_);
 	arrcpy(thearray+i, thearray+i+1, size_-i-1);
@@ -192,7 +197,7 @@ public:
 	set_size(size_ + src.size_);
 	arrcpy(thearray+s,src.thearray, src.size_);	
     }
-    Array<T> subvec(int lower, int upper) {
+    Array<T> slice(int lower, int upper) {
 	assert(lower >= 0 && lower <=upper&& upper <= size_);
 	Array<T> r;
 	int s =upper-lower;
diff --git a/flower/include/vector.hh b/flower/include/vector.hh
index fdeab97f13..7f75735344 100644
--- a/flower/include/vector.hh
+++ b/flower/include/vector.hh
@@ -23,6 +23,11 @@ public:
 	dat.set_size(n);
 	fill(0);
     }
+    void set_dim(int i)
+    {
+	dat.set_size(i);
+    }
+	
     void insert(Real v, int i) {
 	dat.insert(v,i);
     }
diff --git a/flower/include/virtual-methods.hh b/flower/include/virtual-methods.hh
index 88dad44aeb..ceb12694e3 100644
--- a/flower/include/virtual-methods.hh
+++ b/flower/include/virtual-methods.hh
@@ -24,10 +24,18 @@ int a_stupid_nonexistent_function_to_allow_the_semicolon_come_out()
 #define IMPLEMENT_STATIC_NAME(c)\
     char const *c::static_name() { return #c; } 
 
+
+
 #define VIRTUAL_COPY_CONS(T, R)\
   virtual R *clone() const { return  new T(*this); } \
   int  yet_another_stupid_function_to_allow_semicolon()
 
+
+#define DECLARE_VIRTUAL_COPY_CONS(T,R)\
+      virtual R *clone() const
+#define IMPLEMENT_VIRTUAL_COPY_CONS(T,R)\
+	R *T::clone() const { return  new T(*this); } \
+    
 #define IMPLEMENT_IS_TYPE_B(D) 							   \
     IMPLEMENT_STATIC_NAME(D)\
   bool D::static_is_type_b(const char *s)					   \
diff --git a/flower/matrix-storage.cc b/flower/matrix-storage.cc
index 76b89f7f36..54deb7d673 100644
--- a/flower/matrix-storage.cc
+++ b/flower/matrix-storage.cc
@@ -19,8 +19,8 @@ Matrix_storage::set_addition_result(Matrix_storage *&dat, Matrix_storage *right)
 
 	if ( R->band_size_i() > L->band_size_i()) { 
 	    L->set_band_size(R->band_size_i());
-	    return ;
 	}
+	    return ;
     }
     if  (!dat || !dat->is_type_b(Full_storage::static_name() )) {
 
diff --git a/flower/matrix.cc b/flower/matrix.cc
index 78e3d777b0..7e1854827b 100644
--- a/flower/matrix.cc
+++ b/flower/matrix.cc
@@ -108,6 +108,10 @@ Matrix::operator=(Matrix const &m)
 int
 Matrix::band_i()const
 {
+    if ( band_b() ) {
+	Diagonal_storage const * diag = (Diagonal_storage*) dat;
+	return diag->band_size_i();
+    }
     int starty = dim();
     while (starty >= 0 ) {
 	for ( int i = starty, j = 0; i < dim(); i++, j++ )
@@ -206,9 +210,13 @@ operator /(Matrix const& m1,Real a)
     return  m;
 }
 
+/*
+  ugh. Only works for square matrices.
+ */
 void
 Matrix::transpose()		// delegate to storage?
 {
+#if 1
     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
 	if (i >= j)
 	    continue;
@@ -216,6 +224,7 @@ Matrix::transpose()		// delegate to storage?
 	dat->elem(i,j) = dat->elem(j,i);
 	dat->elem(j,i)=r;
     }
+#endif
 }
 
 Matrix
@@ -252,6 +261,7 @@ Matrix::set_product(Matrix const &m1, Matrix const &m2)
     
     if (m1.dat->try_right_multiply(dat, m2.dat))
 	return; 
+    
     for (int i=0, j=0; dat->mult_ok(i,j);
 	 dat->mult_next(i,j)) {
 	Real r=0.0;
diff --git a/flower/rational.cc b/flower/rational.cc
new file mode 100644
index 0000000000..b1056c2430
--- /dev/null
+++ b/flower/rational.cc
@@ -0,0 +1,18 @@
+/*
+  rational.cc -- implement Rational related functions
+
+  source file of the Flower Library
+
+  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+#include "rational.hh"
+#include "string.hh"
+
+void
+print_rat(Rational const &m)
+{
+    cout << String(m) << flush;
+}
+    
+
diff --git a/flower/test/mat-test.cc b/flower/test/mat-test.cc
index f7ce01bfc5..8fdc7d794e 100644
--- a/flower/test/mat-test.cc
+++ b/flower/test/mat-test.cc
@@ -47,7 +47,7 @@ matrix()
     Matrix hilbert(N,N), h2(hilbert);
     for (int i=0; i < N; i++) {
 	for (int j=0; j < N; j++) {
-	    hilbert(i,j) = 1/(i+j+1);
+	    hilbert(i,j) = 1/Real(i+j+1);
 	     h2 (i,j) = (abs(i-j) > 3) ?0 : hilbert(i,j);
 	}
     }
@@ -55,6 +55,10 @@ matrix()
     Choleski_decomposition ch(h2);
     cout << "red Hilbert  " <<  h2;
     cout << "choleski " << ch.L;
+    Matrix T =ch.L.transposed();
+    cout << "L^T " <<  T;
+    cout << "L * L^T" << ch.L * T;
+    cout << "H2^{-1} * H2" << h2 * ch.inverse();
 }
 
 ADD_TEST(matrix);