release: 0.0.2

28 files changed, 1838 insertions, 42 deletions

diff --git a/flower/Makefile b/flower/Makefile
index 5fb4f2a5b6..9db4fbf5e8 100644
--- a/flower/Makefile
+++ b/flower/Makefile
@@ -1,37 +1,31 @@
 MAJVER=1
 MINVER=0
-PATCHLEVEL=1
+PATCHLEVEL=2
 
 PACKAGENAME=flower
 VERSION=$(MAJVER).$(MINVER).$(PATCHLEVEL)
 DNAME=$(PACKAGENAME)-$(VERSION)
-CXXFLAGS+=-g -Wall
-
-cc=lgetopt.cc 	 string.cc dataf.cc textdb.cc unionfind.cc 
-
-templatecc=cursor.cc list.cc 
-inl=findcurs.inl link.inl list.inl
-hh=cursor.hh cursor.inl lgetopt.hh link.hh list.hh \
-	string.hh stringutil.hh vray.hh textdb.hh textstr.hh  assoc.hh\
-	findcurs.hh unionfind.hh compare.hh handle.hh
-
+#DEFINES=-DNDEBUG
+CXXFLAGS+=$(DEFINES) -g -Wall -W -pedantic
 
+include Sources.make
 
 obs=$(cc:.cc=.o)
 staticlib=libflower.a
 $(staticlib): $(obs)
-	ar cr libflower.a $(obs)
+	$(AR) cr libflower.a $(obs)
 
 include depend
 
-depend: Makefile
+depend: Sources.make
 	$(CXX) -MM $(cc) > depend
 
 
 clean:
-	rm depend $(obs) $(staticlib)
-
-DFILES=$(hh) $(cc) $(inl) $(templatecc) Makefile
+	rm -f $(obs) $(staticlib)
+realclean: clean
+	rm -f depend
+DFILES=$(hh) $(cc) $(inl) $(templatecc) Makefile Sources.make TODO
 DDIR=$(DNAME)
 
 dist:
@@ -39,5 +33,6 @@ dist:
 	ln $(DFILES) $(DDIR)/
 	tar cfz $(DNAME).tar.gz $(DDIR)/*
 	rm -rf $(DDIR)/
+TAGS:
+	etags -CT $(inl) $(cc) $(hh)
 
-	
diff --git a/flower/Sources.make b/flower/Sources.make
new file mode 100644
index 0000000000..18fb18cb24
--- /dev/null
+++ b/flower/Sources.make
@@ -0,0 +1,12 @@
+
+cc=lgetopt.cc 	 string.cc dataf.cc textdb.cc unionfind.cc  \
+	smat.cc matrix.cc choleski.cc vector.cc dstream.cc\
+	matdebug.cc
+
+templatecc=cursor.cc list.cc tsmat.cc 
+inl=findcurs.inl link.inl list.inl
+hh=cursor.hh cursor.inl lgetopt.hh link.hh list.hh dstream.hh \
+	string.hh stringutil.hh vray.hh textdb.hh textstr.hh  assoc.hh\
+	findcurs.hh unionfind.hh compare.hh handle.hh matrix.hh\
+	smat.hh vsmat.hh  vector.hh  real.hh choleski.hh\
+	tsmat.hh tvsmat.hh
diff --git a/flower/TODO b/flower/TODO
new file mode 100644
index 0000000000..c8f3e5ad52
--- /dev/null
+++ b/flower/TODO
@@ -0,0 +1,20 @@
+
+	* change String::pos
+
+		s[s.pos('%')] == '%'
+
+	would be nice
+
+	* use template handle in handle.hh for strings.
+
+	* Restricted cursor/list: make sublist from a list, and use rcursor as if list is as big as the sublist.
+
+	* Cursor signedcompare
+
+	* int Cursor::op-(Cursor)
+
+	* move towards gnu?
+
+		parsestream.h
+		vector.h
+
diff --git a/flower/assoc.hh b/flower/assoc.hh
index a2173026f4..84a54c9a72 100644
--- a/flower/assoc.hh
+++ b/flower/assoc.hh
@@ -70,4 +70,6 @@ public:
     }
 
 };
+/** mindblowingly stupid Associative array implementation
+ */
 #endif
diff --git a/flower/choleski.cc b/flower/choleski.cc
new file mode 100644
index 0000000000..82740a27c1
--- /dev/null
+++ b/flower/choleski.cc
@@ -0,0 +1,97 @@
+#include "choleski.hh"
+const Real EPS = 1e-7;		// so sue me. Hard coded
+
+Vector
+Choleski_decomposition::solve(Vector rhs)const
+{
+    int n= rhs.dim();
+    assert(n == L.dim());
+    Vector y(n);
+
+    // forward substitution
+    for (int i=0; i < n; i++) {
+	Real sum(0.0);
+	for (int j=0; j < i; j++)
+	    sum += y(j) * L(i,j);
+	y(i) = (rhs(i) - sum)/L(i,i);
+    }
+    for (int i=0; i < n; i++) {
+	assert(D(i));
+	y(i) /= D(i);
+    }
+
+    // backward subst
+    Vector x(n);
+    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;
+}
+
+/*
+  Standard matrix algorithm.
+  */
+
+Choleski_decomposition::Choleski_decomposition(Matrix P)
+    : L(P.dim()), D(P.dim())
+{
+    int n = P.dim();
+    assert((P-P.transposed()).norm()/P.norm() < EPS);
+
+    L.unit();
+    for (int k= 0; k < n; k++) {
+	for (int j = 0; j < k; j++){
+	    Real sum(0.0);
+	    for (int l=0; l < j; l++)
+		sum += L(k,l)*L(j,l)*D(l);
+	    L(k,j) = (P(k,j) - sum)/D(j);
+	}
+	Real sum=0.0;
+	
+	for (int l=0; l < k; l++)
+	    sum += sqr(L(k,l))*D(l);
+	Real d = P(k,k) - sum;
+	D(k) = d;
+    }
+
+#ifdef NDEBUG
+    assert((original()-P).norm() / P.norm() < EPS);
+#endif
+}
+     
+Matrix
+Choleski_decomposition::original() const
+{
+    Matrix T(L.dim());
+    T.set_diag(D);
+    return L*T*L.transposed();
+}
+
+Matrix
+Choleski_decomposition::inverse() const
+{
+    int n=L.dim();
+    Matrix invm(n);
+    Vector e_i(n);
+    for (int i = 0; i < n; i++) {
+	e_i.set_unit(i);
+	Vector inv(solve(e_i));
+	for (int j = 0 ; j<n; j++)
+	    invm(i,j) = inv(j);
+    }
+    
+#ifdef NDEBUG
+    Matrix I1(n), I2(original());
+    I1.unit();
+    assert((I1-original()*invm).norm()/original.norm() < EPS);
+#endif
+    
+    return invm;
+}
+
+
+
+
diff --git a/flower/choleski.hh b/flower/choleski.hh
new file mode 100644
index 0000000000..c6cb91723b
--- /dev/null
+++ b/flower/choleski.hh
@@ -0,0 +1,46 @@
+#ifndef CHOLESKI_HH
+#define CHOLESKI_HH
+
+#include "matrix.hh"
+
+struct Choleski_decomposition {
+
+    /// lower triangle of Choleski decomposition
+    Matrix L;
+
+    /// diagonal 
+    Vector D;
+    ///Create decomposition of P
+    Choleski_decomposition(Matrix P);
+    /**
+    PRE
+    P needs to be symmetric positive definite
+    */
+    
+    Vector solve(Vector rhs) const;
+    Vector operator * (Vector rhs) const { return solve (rhs); }
+    /**
+    solve Px = rhs
+    */
+
+    Matrix inverse() const;
+    /**
+    return the inverse of the matrix P.
+    */
+
+    Matrix original() const;
+    /**
+    return P,  calc'ed from L and D
+    */
+	 
+};
+/**
+    structure for using the LU decomposition of a positive definite .
+
+    #P# is split  into
+
+    LD transpose(L)
+    */
+    
+
+#endif
diff --git a/flower/cursor.hh b/flower/cursor.hh
index b053591dd4..45c02aea15 100644
--- a/flower/cursor.hh
+++ b/flower/cursor.hh
@@ -70,9 +70,13 @@ class Cursor
       cursor points to same object, cursor.previous()
       is newly inserted object.
       */
-    /// remove and cleanup Link // HWN: backspace or del?
-    void remove();		
+    
+    ///
+    void backspace();
 
+    /// 
+    void del();
+    
     /// access the list this came from
     const List<T>& list() const ;
     Link<T>* pointer();
diff --git a/flower/cursor.inl b/flower/cursor.inl
index 0243d14a2b..a4e03492b5 100644
--- a/flower/cursor.inl
+++ b/flower/cursor.inl
@@ -1,4 +1,4 @@
- // cursor.inl
+ // cursor.inl -*-c++-*-
 #ifndef CURSOR_INL
 #define CURSOR_INL
 #include <assert.h>
@@ -10,7 +10,8 @@ Cursor<T>::Cursor( List<T>& list, Link<T>* pointer ) :
     list_( list )
 {
     if ( list.size() )
-        pointer_ = pointer ? pointer : list.top().pointer_;
+        pointer_ = pointer ? pointer : list.top_;
+    //list.top().pointer_; // ARGH! recursion.
     else
         pointer_ = pointer;
 }
@@ -56,13 +57,24 @@ Cursor<T>::insert( const T& thing )
 
 template<class T>
 inline void
-Cursor<T>::remove()
+Cursor<T>::backspace()
 {
-    assert( pointer_ );
+    Cursor<T> c(*this);
+    c--;        
     list_.remove( *this );
 }
 
 template<class T>
+inline void
+Cursor<T>::del()
+{
+    Cursor<T> c(*this);
+    c++;
+    list_.remove( *this );    
+    *this = c;
+}
+
+template<class T>
 inline const List<T>&
 Cursor<T>::list() const
 {
@@ -97,4 +109,4 @@ Cursor<T>::ok()
     return ( pointer_ != 0 );
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/flower/dstream.cc b/flower/dstream.cc
new file mode 100644
index 0000000000..d8cff69041
--- /dev/null
+++ b/flower/dstream.cc
@@ -0,0 +1,123 @@
+#include <fstream.h>
+
+#include "dstream.hh"
+#include "string.hh"
+#include "textdb.hh"
+
+/// indent of each level 
+const INDTAB = 3;
+
+/*
+  should use Regexp library.
+  */
+static String
+strip_pretty(String pret)
+{
+    String cl(pret.left(pret.pos('(')-1));
+    int l = cl.lastPos(' ');
+    cl = cl.right(cl.len() -l);
+    return cl;
+}
+
+static String
+strip_member(String pret)
+{
+    String cl(pret.left(pret.lastPos(':')-2));
+    return cl;
+}
+
+Dstream&
+Dstream::identify_as(String name)
+{
+    if (!os)
+	return *this;
+    
+    String mem(strip_pretty(name));
+    String cl(strip_member(mem));
+    String idx = cl;
+    
+    if (silent.elt_query(mem))
+	idx  = mem;
+    else if (silent.elt_query(cl))
+	idx = cl;
+    else {
+	silent[idx] = false;
+    }
+    local_silence = silent[idx];
+    if (classname != idx && !local_silence) {
+	classname=idx;
+	*os << "[" << classname << ":]";
+    }
+    return *this;
+}
+
+bool
+Dstream::silence(String s)
+{
+    if (!silent.elt_query(s))
+	return false;
+    return silent[s];
+}
+///
+Dstream &
+Dstream::operator<<(String s)
+{
+    if (local_silence|| !os)
+	return *this;
+    
+    for (const char *cp = s  ; *cp; cp++)
+	switch(*cp) 
+	    {
+	    case '{':
+	    case '[':
+	    case '(': indentlvl += INDTAB;
+		*os << *cp;		
+		break;
+		
+	    case ')':
+	    case ']':
+	    case '}':
+		indentlvl -= INDTAB;
+		*os << *cp		;
+		
+		assert  (indentlvl>=0) ;
+		break;
+		
+	    case '\n':
+		*os << '\n' << String (' ', indentlvl) << flush;
+		break;	      
+	    default:
+		*os << *cp;
+		break;
+	    }
+    return *this;    
+}
+
+/** only output possibility. Delegates all conversion to String class.
+ */
+
+Dstream::Dstream(ostream *r, const char * cfg_nm )
+{
+    os = r;
+    if (!os)
+	return;
+    indentlvl = 0;
+    
+    const char * fn =cfg_nm ? cfg_nm : ".dstreamrc";
+    {
+	ifstream ifs(fn);	// can't open
+	if (!ifs)
+	    return;
+    }
+    //    cerr << "(" << fn;
+    Text_db cfg(fn);
+    while (! cfg.eof()){	     
+	 Text_record  r(  cfg++);
+	 assert(r.sz() == 2);
+	 silent[r[0]] = r[1].to_bool();
+    }
+    //  cerr <<")";
+}
+
+
+
diff --git a/flower/dstream.hh b/flower/dstream.hh
new file mode 100644
index 0000000000..72d0897201
--- /dev/null
+++ b/flower/dstream.hh
@@ -0,0 +1,48 @@
+// debug_stream
+
+#ifndef DSTREAM_HH
+#define DSTREAM_HH
+
+#include "string.hh"
+#include "assoc.hh"
+
+const char eol= '\n';
+
+/// debug stream
+class Dstream
+{
+    ostream *os;
+    int indentlvl;
+    bool local_silence;
+    String classname;
+
+    Assoc<String, bool> silent;
+public:
+
+    bool silence(String);
+    
+    Dstream(ostream *r, const char  * rcfile);
+    /**
+      if rcfile == 0, then do not read any rc file
+      */
+      
+    Dstream &identify_as(String s);
+    
+    Dstream &operator << (String s);
+};
+ /**
+   a class for providing debug output of nested structures,
+   with indents according to \{\}()[].
+
+   One can turn on and off specific messages using the Assoc silent.
+   This can be done automatically:
+
+   #define DEBUG  dstream_.identify_as(__PRETTY_FUNCTION__) 
+
+   DEBUG << "a message\n";
+   
+   Init for the class names which should be silent can be given in a rc file. 
+   
+  */
+#endif
+
diff --git a/flower/link.inl b/flower/link.inl
index e76caab4f1..3926d6bc2a 100644
--- a/flower/link.inl
+++ b/flower/link.inl
@@ -7,12 +7,14 @@ inline
 void
 Link<T>::OK() const
 {
+#ifndef NDEBUG
     if (previous_) {
 	assert(previous_->next_ == this);
     }
     if (next_) {
 	assert(next_->previous_ == this);
     }
+#endif    
 }
 
 template<class T>
diff --git a/flower/list.cc b/flower/list.cc
index f93352f507..782b3da86e 100644
--- a/flower/list.cc
+++ b/flower/list.cc
@@ -37,6 +37,7 @@ List<T>::top()
 	assert( t != top_ );	// silly link
 	while ( t )
 	    {
+		assert(false);	// this is even more silly.
 	    top_ = t;
 	    t = top_->previous();
 	    }
diff --git a/flower/list.hh b/flower/list.hh
index 97e672f2e0..a45252b8c3 100644
--- a/flower/list.hh
+++ b/flower/list.hh
@@ -35,7 +35,12 @@ class List
     /// put thing before #before_me#
     void insert( const T& thing, Cursor<T> before_me );
     virtual void remove( Cursor<T> me );
+    /**
+      Remove link pointed to by me.
 
+      POST
+      none; WARNING: do not use #me#.
+     */
     int size_;
     Link<T>* top_;
     Link<T>* bottom_;
diff --git a/flower/list.inl b/flower/list.inl
index 54df457481..8396156b6a 100644
--- a/flower/list.inl
+++ b/flower/list.inl
@@ -22,8 +22,12 @@ template<class T>
 inline
 List<T>::~List()
 {
-    for ( Cursor<T> c( *this ); c.forward(); c++ )
-        remove( c );
+    Cursor<T> next(*this);
+    for ( Cursor<T> c( *this ); c.ok(); c = next ) {
+	next = c;
+	next++;
+	remove( c );
+    }
 }
 
 template<class T>
@@ -94,12 +98,12 @@ template<class T>
 inline void
 List<T>::remove( Cursor<T> me )
 {
-    if ( me.ok() )
-	{
-	me.pointer()->remove(*this);
-	delete me.pointer();
+    if ( me.ok() ){
+	Link<T> *lp = me.pointer();	
+	lp->remove(*this);
+	delete lp;
         size_--;
-	}
+    }
 }
 
 template<class T>
@@ -127,8 +131,12 @@ template<class T>
 inline
 PointerList<T>::~PointerList()
 {
-    for ( Cursor<T> c( *this ); c.forward(); c++ )
-        remove( c );
+    Cursor<T> next(*this);
+    for ( Cursor<T> c( *this ); c.ok(); c = next ) {
+	next = c;
+	next++;
+	remove( c );		// PointerList::remove deletes the  real data
+    }
 }
 
 template<class T>
@@ -136,7 +144,7 @@ inline void
 PointerList_print( PointerList<T> const & l  ) 
 {
     List<T>& promises_to_be_const = (List<T>&) l;
-    for ( Cursor<T> c( promises_to_be_const ); c.forward(); c++ )
+    for ( Cursor<T> c( promises_to_be_const ); c.ok(); c++ )
         (*c)->print();  
 }
 
diff --git a/flower/matdebug.cc b/flower/matdebug.cc
new file mode 100644
index 0000000000..5a31312386
--- /dev/null
+++ b/flower/matdebug.cc
@@ -0,0 +1,55 @@
+#include "dstream.hh"
+#include "matrix.hh"
+
+static Dstream *dout = new Dstream(0,0);
+
+void set_matrix_debug(Dstream&ds)
+{
+    dout = &ds;
+}
+
+Matrix::operator String() const
+{
+    String s("matrix {\n");
+#ifndef NPRINT
+    for (int i=0; i< rows(); i++){
+	for (int j = 0; j < cols(); j++) {
+	    s+= String(dat->elem(i,j), "%6f ");
+	}
+	s+="\n";
+    }
+    s+="}\n";
+#endif
+    return s;
+}
+
+
+void
+Matrix::print() const
+{
+#ifndef NPRINT
+    *dout << *this;
+#endif
+}
+
+Vector::operator String() const
+{
+    int i=0;
+    String s("vector [");
+#ifndef NDEBUG
+    for (; i < dim(); i++) {
+	s += String(dat[i], "%6f") + ' ';
+    }
+#endif
+    s+="]";
+    return s;
+}
+
+
+void
+Vector::print() const
+{
+#ifndef NDEBUG
+    *dout << *this<<'\n';
+#endif
+}
diff --git a/flower/matrix.cc b/flower/matrix.cc
new file mode 100644
index 0000000000..1afe5951c8
--- /dev/null
+++ b/flower/matrix.cc
@@ -0,0 +1,260 @@
+#include "matrix.hh"
+#include "string.hh"
+
+
+Real
+Matrix::norm() const
+{
+    Real r =0.0;
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+  	r += sqr(dat->elem(i,j));
+    return sqrt(r);
+}
+
+//inline
+Real
+Matrix::operator()(int i,int j) const
+{
+    assert(i >= 0 && j >= 0);
+    assert(i < rows() && j < cols());
+    return dat->elem(i,j);
+}
+
+//inline
+Real &
+Matrix::operator()(int i, int j)
+{
+    assert(i >= 0 && j >= 0);
+    assert(i < rows() && j < cols());
+    return dat->elem(i,j);
+}
+
+void
+Matrix::fill(Real r)
+{
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+	dat->elem(i,j)=r;
+}
+
+void
+Matrix::set_diag(Real r)
+{
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))	  
+ 	dat->elem(i,j)=(i==j) ? r: 0.0;
+}
+
+void
+Matrix::set_diag(Vector d)
+{
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))	  
+ 	dat->elem(i,j)=(i==j) ? d(i): 0.0;
+}
+
+void
+Matrix::operator+=(const Matrix&m)
+{
+    assert(m.cols() == cols());
+    assert(m.rows() == rows());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+  	dat->elem(i,j) += m(i,j);
+}
+ 
+void
+Matrix::operator-=(const Matrix&m)
+{
+    assert(m.cols() == cols());
+    assert(m.rows() == rows());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+  	dat->elem(i,j) -= m(i,j);
+}
+
+
+void
+Matrix::operator*=(Real a)
+{
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+	dat->elem(i,j) *= a;
+}
+
+void
+Matrix::operator=(const Matrix&m)
+{
+    if (&m == this)
+	return ;
+    delete dat;
+    dat = m.dat->clone();
+}
+    
+Matrix::Matrix(const Matrix &m)
+{
+    m.OK();
+    
+    dat = m.dat->clone();
+}
+
+
+Matrix::Matrix(int n, int m)
+{
+    dat = virtual_smat::get_full(n,m);
+    fill(0);
+}
+
+Matrix::Matrix(int n)
+{
+    dat = virtual_smat::get_full(n,n);
+    fill(0);
+}
+
+Matrix::Matrix(Vector v, Vector w)
+{   
+    dat = virtual_smat::get_full(v.dim(), w.dim());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+	dat->elem(i,j)=v(i)*w(j);
+}
+
+
+Vector
+Matrix::row(int k) const
+{
+    int n=cols();
+
+    
+    Vector v(n);
+    for(int i=0; i < n; i++)
+	v(i)=dat->elem(k,i);
+
+    return v;
+}
+
+Vector
+Matrix::col(int k) const
+{
+    int n=rows();
+    Vector v(n);
+    for(int i=0; i < n; i++)
+	v(i)=dat->elem(i,k);
+    return v;
+}
+
+Vector
+Matrix::left_multiply(const Vector& v) const
+{
+     Vector dest(v.dim());
+    assert(dat->cols()==v.dim());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+	dest(i)+= dat->elem(j,i)*v(j);
+    return dest;
+}
+
+Vector
+Matrix::operator *(const Vector& v) const
+{
+    Vector dest(rows());
+    assert(dat->cols()==v.dim());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j))
+	dest(i)+= dat->elem(i,j)*v(j);
+    return dest;
+}
+
+Matrix
+operator /(Matrix const& m1,Real a)
+{
+    Matrix m(m1);
+    m /= a;
+    return  m;
+}
+
+void
+Matrix::transpose()		// delegate to storage?
+{
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
+	if (i >= j)
+	    continue;
+	Real r=dat->elem(i,j);
+	dat->elem(i,j) = dat->elem(j,i);
+	dat->elem(j,i)=r;
+    }
+}
+
+Matrix
+Matrix::operator-() const
+{
+    OK();
+    Matrix m(*this);
+    m*=-1.0;    
+    return m;
+}
+
+Matrix
+Matrix::transposed() const
+{
+    Matrix m(*this);
+    m.transpose();
+    return m;
+}
+
+
+/* should do something smarter: bandmatrix * bandmatrix is also banded matrix.  */
+Matrix
+operator *(const Matrix &m1, const Matrix &m2)
+{
+    Matrix result(m1.rows(), m2.cols());
+    result.set_product(m1,m2);
+    return result;
+}
+
+void
+Matrix::set_product(const Matrix &m1, const Matrix &m2)
+{
+    assert(m1.cols()==m2.rows());
+    assert(cols()==m2.cols() && rows()==m1.rows());
+    
+    for (int i=0, j=0; dat->mult_ok(i,j);
+	 dat->mult_next(i,j)) {
+	Real r=0.0;
+	for (int k = 0; k < m1.cols(); k++)
+	    r += m1(i,k)*m2(k,j);
+	dat->elem(i,j)=r;
+    }
+}
+
+void
+Matrix::insert_row(Vector v, int k)
+{
+    assert(v.dim()==cols());
+    dat->insert_row(k);
+    for (int j=0; j < cols(); j++)
+	dat->elem(k,j)=v(j);
+}
+
+
+void
+Matrix::swap_columns(int c1, int c2)
+{
+    assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
+    for (int i=0; i< rows(); i++) {
+	Real r=dat->elem(i,c1);
+	dat->elem(i,c1) = dat->elem(i,c2);
+	dat->elem(i,c2)=r;
+    }
+}
+
+void
+Matrix::swap_rows(int c1, int c2)
+{
+    assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0);
+    for (int i=0; i< cols(); i++) {
+	Real r=dat->elem(c1,i);
+	dat->elem(c1,i) = dat->elem(c2,i);
+	dat->elem(c2,i)=r;
+    }
+}
+
+
+int
+Matrix::dim() const
+{
+    assert(cols() == rows());
+    return rows();
+}
+
diff --git a/flower/matrix.hh b/flower/matrix.hh
new file mode 100644
index 0000000000..e092a5fbe6
--- /dev/null
+++ b/flower/matrix.hh
@@ -0,0 +1,141 @@
+#ifndef MATRIX_HH
+#define MATRIX_HH
+
+
+#include "vsmat.hh"
+#include "vector.hh"
+
+/// a Real matrix
+class Matrix {
+    virtual_smat *dat;
+    
+public:
+    void OK() const { dat->OK(); }
+    int cols() const { return dat->cols(); }
+    int rows() const { return dat->rows(); }
+
+    /// return the size of a matrix
+    int dim() const;
+    /**
+      PRE
+      the matrix needs to be square.
+    */
+     
+    // Matrix() { dat = 0; } 
+    ~Matrix() { delete dat; }
+
+    /// set entries to r 
+    void fill(Real r);
+
+    /// set diagonal to d
+    void set_diag(Real d);
+
+    void set_diag(Vector d);
+    /// set unit matrix
+    void unit() { set_diag(1.0); }
+
+    void operator+=(const Matrix&m);
+        void operator-=(const Matrix&m);    
+    void operator*=(Real a);
+    void operator/=(Real a) { (*this) *= 1/a; }
+    
+    /// add a row
+    void insert_row(Vector v,int k);
+    /**
+      add a row to the matrix before  row k
+
+      PRE
+      v.dim() == cols()
+      0 <= k <= rows()
+    */
+    ///
+    void delete_row(int k) { dat->delete_row(k); }
+    /**
+      delete a row from this matrix.
+
+      PRE
+      0 <= k < rows();
+    */
+    void delete_column(int k) { dat->delete_column(k); }
+    ///
+    Matrix(int n);
+    /**
+      square n matrix, initialised to null
+    */
+    ///
+    Matrix(int n, int m);
+    /**
+      n x m matrix, init to 0
+    */
+    Matrix(const Matrix &m);
+
+    /// dyadic product: v * w.transpose
+    Matrix(Vector v, Vector w);
+    void operator=(const Matrix&m);
+
+    /// access an element
+    Real operator()(int i,int j) const;
+
+    /// access an element
+    Real &operator()(int i, int j);
+
+    /// Matrix multiply with vec (from right)
+    Vector operator *(const Vector &v) const;
+
+    /// set this to m1*m2.
+    void set_product(const Matrix &m1, const Matrix &m2);
+
+
+    Vector left_multiply(Vector const &) const;
+    
+    Matrix operator-() const;
+    
+    /// transpose this.
+    void transpose();
+    
+    /// return a transposed copy.
+    Matrix transposed() const ;
+
+    Real norm() const;
+    /// swap
+    void swap_columns(int c1, int c2);
+    /**
+      PRE
+      0 <= c1,c2 < cols()
+    */
+
+    /// swap
+    void swap_rows(int c1, int c2);
+    /**
+      PRE
+      0 <= c1,c2 < rows()
+    */
+
+
+    Vector row(int ) const;
+    Vector col(int) const;
+
+    operator String() const;
+    void print() const;
+};
+
+/** This is a class for a nonsquare block of #Real#s.  The
+    implementation of sparse matrices is done in the appropriate #smat#
+    class. Matrix only does the mathematical actions (adding,
+    multiplying, etc.)
+
+    
+    TODO
+    implement ref counting?  */
+
+
+inline Vector
+operator *(Vector &v, const Matrix& m) { return m.left_multiply(v); }
+Matrix operator *(const Matrix& m1,const Matrix &m2);
+Matrix operator /(const Matrix &m1,Real a);
+inline Matrix operator -(Matrix m1,const Matrix m2)
+{
+    m1 -= m2;
+    return m1;
+}
+#endif
diff --git a/flower/real.hh b/flower/real.hh
new file mode 100644
index 0000000000..1f2187c8c1
--- /dev/null
+++ b/flower/real.hh
@@ -0,0 +1,25 @@
+#ifndef REAL_HH
+#define REAL_HH
+typedef double Real;
+inline Real sqr(Real x){
+    return x*x;
+}
+inline Real MIN(Real x, Real y) {
+    return (x < y)? x : y;
+}
+
+inline Real MAX(Real x, Real y) {
+    return (x > y) ? x : y;
+}
+
+inline Real ABS(Real x)
+{
+    return (x>0)? x:-x;
+}
+inline
+int sgn(Real x) {
+    if (!x)return 0;
+    return (x > 0) ?1: -1;
+}
+
+#endif
diff --git a/flower/smat.cc b/flower/smat.cc
new file mode 100644
index 0000000000..3830fd6825
--- /dev/null
+++ b/flower/smat.cc
@@ -0,0 +1,187 @@
+#include "smat.hh"
+
+void
+Full_storage::operator=(Full_storage const &fs)
+{
+    resize(fs.h, fs.w);
+    OK();
+    fs.OK();
+    for (int i=0; i<h; i++)
+	for (int j=0; j<w; j++)
+	    els[i][j]= fs.els[i][j];
+}
+
+void
+Full_storage::OK() const
+{
+    //    static Real dummy;		
+    assert(maxh >= h && maxw >= w);
+    assert(h >= 0 && w >= 0);
+    assert(els||!maxh);
+    if (maxh>0) {		// access outer elts.
+	Real *r = els[maxh -1];
+	if (maxw>0) {
+	    assert(r);
+	    Real s = r[maxw -1];
+	    s = sin(s);
+	}
+    }
+}
+void
+Full_storage::resize_cols(int newh)
+{
+    if (newh <= maxh) {
+	h=newh;
+	return;
+    }
+   
+    Real ** newa=new Real*[newh];
+    int j=0;
+    for (; j < h; j++)
+	newa[j] = els[j];
+    for (; j < newh; j++)
+	newa[j] = new Real[maxw];
+    delete[] els;
+    els=newa;
+
+    h = maxh = newh;
+}
+
+void
+Full_storage::resize_rows(int neww)
+{
+    if (neww <= maxw) {
+	w=neww;
+	return;
+    }
+    for (int i=0; i < maxh ; i++) {
+	Real* newa=new Real[neww];
+	for (int k=0; k < w; k++)
+	    newa[k] = els[i][k];
+
+	delete[] els[i];
+	els[i] = newa;
+    }
+    w = maxw = neww;	
+}
+
+Full_storage::~Full_storage() {
+    for (int i=0; i < maxh; i++)
+	delete [] els[i];
+    delete[] els;
+}
+
+void
+Full_storage::resize(int rows, int cols)
+{
+    OK();
+    resize_cols(rows);
+    resize_rows(cols);
+
+}
+
+
+bool
+Full_storage::mult_ok(int i, int j) const
+{
+    return valid(i,j);
+}
+
+bool
+Full_storage::trans_ok(int i, int j) const
+{
+       return valid(i,j);
+} 
+
+
+void
+Full_storage::trans_next(int &i, int &j) const
+{
+    assert(trans_ok(i,j));
+    i++;
+    if (i >= h) {
+	i=0;
+	j ++;
+    }
+}
+
+void
+Full_storage::mult_next(int &i, int &j) const
+{
+    assert(mult_ok(i,j));
+    j++;
+    if (j >= w) {
+	j=0;
+	i++;
+    }
+}
+
+void
+Full_storage::delete_column(int k)
+{
+    assert(0 <= k &&k<w);    
+    for (int i=0; i< h ; i++)
+	for (int j=k+1; j <w; j++)
+	    els[i][j-1]=els[i][j];
+    w--;
+}
+void
+Full_storage::delete_row(int k)
+{
+    assert(0 <= k &&k<h);
+    for (int i=k+1; i < h ; i++)
+	for (int j=0; j < w; j++)
+	    els[i-1][j]=els[i][j];
+    h--;
+}
+
+
+void
+Full_storage::insert_row(int k)
+{
+    assert(0 <= k&& k <=h);
+    resize_cols(h+1);
+    for (int i=h-1; i > k ; i--)
+	for (int j=0; j <w; j++)
+	    els[i][j]=els[i-1][j];
+
+}
+
+
+svec<Real>
+Full_storage::row(int n) const
+{
+    svec<Real> r;
+    for (int j = 0; j < w; j++)
+	r.add(els[n][j]);
+    return r;
+}
+
+svec<Real>
+Full_storage::column(int n) const
+{
+    
+    svec<Real> r;
+    for (int i = 0; i<h; i++)
+	r.add(els[i][n]);
+    return r;
+}
+
+
+Full_storage::Full_storage(Full_storage&s)
+{
+    init();
+    (*this) = s;
+}
+virtual_smat*
+Full_storage::clone()
+{
+    return new Full_storage(*this);
+}
+/****************************************************************/
+
+virtual_smat *
+virtual_smat::get_full(int n, int m)
+{
+    return new Full_storage(n,m);
+}
diff --git a/flower/smat.hh b/flower/smat.hh
new file mode 100644
index 0000000000..c05d032d89
--- /dev/null
+++ b/flower/smat.hh
@@ -0,0 +1,93 @@
+#ifndef SMAT_HH
+#define SMAT_HH
+#include "vray.hh"
+#include "vsmat.hh"
+#include "real.hh"
+/// simplest matrix storage. refer to its baseclass for the doco.
+class Full_storage : public virtual_smat
+{
+    /// height, width
+    int h,w;
+    /// maxima.
+    int maxh, maxw;
+    
+    /// the storage
+    Real** els;
+    void
+    init() {
+	els=0;
+	h=w=maxh=maxw=0;
+
+    }
+    
+    bool valid(int i, int j) const {
+	return (i>=0 && i < h)
+	    && (j < w && j >=0);
+    }
+    
+
+    void resize_rows(int);
+    void resize_cols(int);
+
+public:
+    virtual int rows() const {
+	return h;
+    }
+    virtual int cols() const {
+	return w;
+    }
+    
+    
+    virtual void set_size(int i, int j)
+    {
+	resize(i,j); //this could be more efficient.
+    }
+    
+    virtual void set_size(int i) {
+	set_size(i,i);
+    }
+    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[i][j];
+    }
+    virtual const Real& elem(int i, int j) const {
+	assert(valid(i,j));
+	return els[i][j];
+    }
+    virtual svec<Real> row(int i) const;
+    virtual svec<Real> column(int j) const;
+
+    Full_storage() {
+	init();
+    }
+    Full_storage(int i, int j) {
+	init();
+	set_size(i,j);
+    }
+    Full_storage(Full_storage&);
+    Full_storage(int i) {
+	init();
+	set_size(i);
+    }
+    void OK() const;
+    void operator=(Full_storage const &);
+    
+    virtual void insert_row(int k);
+    virtual void delete_row(int k);
+    virtual void delete_column(int k);
+
+    
+    ~Full_storage();
+    virtual bool mult_ok(int i, int j)const;
+    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 virtual_smat * clone();
+};
+
+#endif
diff --git a/flower/string.cc b/flower/string.cc
index da4941da7b..2e8e8bf5c6 100644
--- a/flower/string.cc
+++ b/flower/string.cc
@@ -15,7 +15,7 @@
 //#include "globals.hh"
 #include "string.hh"
 
-char* strlwr( char* s )
+static char* strlwr( char* s )
 {
     char* p = s;
 
@@ -27,7 +27,7 @@ char* strlwr( char* s )
     return s;
 }
 
-char* strupr( char* s )
+static char* strupr( char* s )
 {
     char* p = s;
 
@@ -320,7 +320,13 @@ String String::lower()
 
 String::String (double f, const char *fmt)
 {
-    char buf[100]; // ugly
+    /* worst case would be printing HUGE (or 1/HUGE), which is approx
+      2e318, this number would have approx 318 zero's in its string.
+
+      1024 is a safe length for the buffer
+      */
+
+    char buf[1024]; 
     if (!fmt)
 	sprintf(buf, "%f", f);
     else
diff --git a/flower/stringutil.hh b/flower/stringutil.hh
index b7503b7d2a..fe0b8e307b 100644
--- a/flower/stringutil.hh
+++ b/flower/stringutil.hh
@@ -2,17 +2,15 @@
 #define STRINGUTIL_HH
 #include <assert.h>
 
-#ifndef EVELYN
+#if !defined(NDEBUG)
 #define NDEBUG BLONDE
-  // switching into "blonde" mode
-  // i trust stringclass nowadays.
 #endif
 
 class String_handle;
 /// Internal String struct
 class StringData {
     // GNU malloc: storage overhead is 8 bytes anyway.
-    const int INITIALMAX = 8;  
+    const int INITIALMAX =8;	// how to do this in ANSI C++ ?
 
 friend class String_handle;
     int maxlen;	// maxlen is arraysize-1
@@ -156,7 +154,8 @@ friend class String_handle;
 
 /**
    the data itself. Handles simple tasks (resizing, resetting)
- */
+   */
+
 /****************************************************************/
 /// ref. counting for strings
 class String_handle {
diff --git a/flower/tsmat.cc b/flower/tsmat.cc
new file mode 100644
index 0000000000..0e5407ae70
--- /dev/null
+++ b/flower/tsmat.cc
@@ -0,0 +1,151 @@
+#include "smat.hh"
+
+template<class T>
+void
+Full_storage<T>::operator=(Full_storage const &fs)
+{
+    resize(fs.h, fs.w);
+    for (int i=0; i<h; i++)
+	for (int j=0; i<w; j++)
+	    els[i][j]= fs.els[i][j];
+}
+
+template<class T>
+void
+Full_storage<T>::OK() const
+{
+    assert(maxh >= h && maxw >= w);
+    assert(h >= 0 && w >= 0);
+}
+template<class T>
+void
+Full_storage<T>::resize_cols(int newh)
+{
+    if (newh <= maxh) {
+	h=newh;
+	return;
+    }
+   
+    T** newa=new T*[newh];
+    int j=0;
+    for (; j < h; j++)
+	newa[j] = els[j];
+    for (; j < newh; j++)
+	newa[j] = new T[w];
+    delete[] els;
+    els=newa;
+    maxh = newh;
+}
+
+template<class T>
+void
+Full_storage<T>::resize_rows(int neww)
+{
+    if (neww <= maxw) {
+	w=neww;
+	return;
+    }
+    for (int i=0; i < h ; i++) {
+	T* newa=new T[neww];
+	for (int k=0; k < w; k++)
+	    newa[k] = els[i][k];
+
+	delete[] els[i];
+	els[i] = newa;
+	maxw = neww;
+    }
+}
+
+template<class T>
+Full_storage<T>::~Full_storage() {
+    for (int i=0; i < maxh; i++)
+	delete [] els[i];
+    delete[] els;
+}
+
+template<class T>
+void
+Full_storage<T>::resize(int i, int j)
+{
+    resize_cols(i);
+    resize_rows(j);
+}
+
+template<class T>
+void
+Full_storage<T>::set_size(int i, int j)
+{    
+    resize(i,j)
+}
+
+template<class T>
+bool
+Full_storage<T>::mult_ok(int i, int j) const
+{
+    return valid(i,j);
+}
+
+template<class T>
+bool
+Full_storage<T>::trans_ok(int i, int j) const
+{
+       return valid(i,j);
+} 
+
+
+template<class T>
+void
+Full_storage<T>::trans_next(int &i, int &j) const
+{
+    assert(trans_ok(i,j));
+    i++;
+    if (i >= h) {
+	i=0;
+	j ++;
+    }
+}
+
+template<class T>
+void
+Full_storage<T>::mult_next(int &i, int &j) const
+{
+    assert(mult_ok(i,j));
+    j++;
+    if (j >= w) {
+	j=0;
+	i++;
+    }
+}
+
+template<class T>
+void
+Full_storage<T>::delete_row(int k)
+{
+    assert(0 <= k <h);
+    for (int i=h-1; i > k ; i++)
+	for (int j=0; j < w; j++)
+	    els[i-1][j]=els[i][j];
+}
+
+template<class T>
+void
+Full_storage<T>::insert_row(int k)
+{
+    assert(0 <= k <=h);
+    resize_cols(h+1);
+    for (int i=h-1; i > k ; i++)
+	for (int j=0; j <w; j++)
+	    els[i][j]=els[i-1][j];
+}
+
+/****************************************************************/
+
+template<class T>
+virtual_smat<T> *
+virtual_smat<T>::get_full(int n, int m)
+{
+    return new Full_storage<T>(n,m);
+}
+#include "real.hh"
+
+template Full_storage<Real>;
diff --git a/flower/tsmat.hh b/flower/tsmat.hh
new file mode 100644
index 0000000000..0fe3f7af79
--- /dev/null
+++ b/flower/tsmat.hh
@@ -0,0 +1,92 @@
+#ifndef SMAT_HH
+#define SMAT_HH
+#include "vray.hh"
+#include "vsmat.hh"
+
+/// simplest matrix storage. refer to its baseclass for the doco.
+template<class T>
+class Full_storage : public virtual_smat<T>
+{
+    /// height, width
+    int h,w;
+    /// maxima.
+    int maxh, maxw;
+    
+    /// the storage
+    T** els;
+    void
+    init() {
+	els=0;
+	h=w=maxh=maxw=0;
+
+    }
+    
+    bool valid(int i, int j) const {
+	return (i>=0 && i < h)
+	    && (j < w && j >=0);
+    }
+    
+
+    void resize_rows(int);
+    void resize_cols(int);
+
+public:
+    virtual int rows() const {
+	return h;
+    }
+    virtual int cols() const {
+	return w;
+    }
+    
+    
+    virtual void set_size(int i, int j)
+    {
+	resize(i,j); //this could be more efficient.
+    }
+    
+    virtual void set_size(int i) {
+	set_size(i,i);
+    }
+    virtual void resize(int i, int j);
+    virtual void resize(int i) {
+	resize(i,i);
+    }
+    
+    virtual T& elem(int i,int j) {
+	assert(valid(i,j));
+	return els[i][j];
+    }
+    virtual const T& elem(int i, int j) const {
+	assert(valid(i,j));
+	return els[i][j];
+    }
+    virtual svec<T> row(int i) const;
+    virtual svec<T> column(int j) const;
+
+    Full_storage() {
+	init();
+    }
+    Full_storage(int i, int j) {
+	init();
+	set_size(i,j);
+    }
+    Full_storage(int i) {
+	init();
+	set_size(i);
+    }
+    void OK() const;
+    void operator=(Full_storage const &);
+    
+    virtual void insert_row(int k);
+    virtual void delete_row(int k);
+
+    
+    ~Full_storage();
+    virtual bool mult_ok(int i, int j)const;
+    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;
+
+};
+
+#endif
diff --git a/flower/tvsmat.hh b/flower/tvsmat.hh
new file mode 100644
index 0000000000..5e79634b55
--- /dev/null
+++ b/flower/tvsmat.hh
@@ -0,0 +1,140 @@
+#ifndef VSMAT_HH
+#define VSMAT_HH
+#include "vray.hh"
+
+/// a  matrix storage baseclass.
+class virtual_smat {
+    
+
+public:
+    /// check invariants
+    virtual void OK() const=0;
+    
+    /// height of matrix
+    virtual int rows() const = 0;
+
+    /// width of matrix
+    virtual int cols() const = 0;
+    
+    /// set the size. contents lost
+    virtual void set_size(int i, int j) = 0;
+    /**      
+      PRE
+      i >=0, j>=0
+    */
+    
+    /// set the size to square dimen. contents lost
+    virtual void set_size(int i) = 0;
+    /**
+      PRE
+      i>=0
+    */
+    /// set the size to i
+    virtual void resize(int i, int j) = 0;
+    /**
+
+      keep contents. If enlarged contents unspecified
+        
+      PRE
+      i>=0, j>=0
+    
+    */
+
+    /// set the size to square dimen. contents kept
+    virtual void resize(int i) = 0;
+    /**    
+    Keep contents. If enlarged contents are unspecified
+    
+    PRE
+    i>=0  
+    */
+    
+    /// access an element
+    virtual Real& elem(int i,int j) = 0;
+    /**
+    access an element.
+
+    Generate an errormessage, if this happens
+    in the 0-part of a sparse matrix.
+    */
+
+    /// access a element, no modify
+    virtual const Real& elem(int i, int j) const = 0;
+
+#if 1
+    virtual svec<Real> row(int i) const = 0;
+    virtual svec<Real> column(int j) const = 0;
+#endif
+
+    /// add a row
+    virtual void insert_row(int k)=0;
+    /**
+    add a row to the matrix before  row k. Contents
+    of added row are unspecified
+
+      0 <= k <= rows()
+    */
+
+    /// delete a row 
+    virtual void delete_row(int k)=0;
+      /**
+      delete a row from this matrix.
+
+      PRE
+      0 <= k < rows();
+    */
+    virtual ~virtual_smat() { }
+    virtual virtual_smat *clone()=0;
+
+
+    /// is there a next?
+    virtual bool mult_ok(int i, int j) const=0;
+    /**
+      at end of matrix? when doing loop
+
+      for(i=0; i<h; i++)
+        for(j=0; j<w; j++)
+          ...
+
+    */
+    /// iterate
+    virtual void mult_next(int &i, int &j) const  = 0;
+    /**
+      walk through matrix (regular multiply)
+      get next j for row i, or get next row i and reset j.
+      this will make sparse matrix implementation easy.
+    
+      PRE
+      mult_ok(i,j)
+     */
+    virtual bool trans_ok(int i, int j) const=0;
+    /**
+      valid matrix entry. return false if at end of row
+    */
+    virtual void trans_next(int &i, int &j) const  = 0;
+    /**
+      walk through matrix (transposed multiply).
+      Get next i (for column j)
+    
+      PRE
+      ver_ok(i,j)
+     */
+
+    /// generate a "Full_storage" matrix    
+    virtual_smat<Real> *get_full(int n, int m);
+
+};
+    
+/** base class for interface with matrix storageclasses.  There are no
+    iterators for matrixclasses, since matrices are (like arrays)
+    explicitly int-indexed.
+
+    Iteration is provided by *_next, *_ok, which update and check both
+    index variables simultaneously.
+
+    TODO
+    determine type of product matrix.
+
+*/
+
+#endif
diff --git a/flower/vector.cc b/flower/vector.cc
new file mode 100644
index 0000000000..2cbfdfc899
--- /dev/null
+++ b/flower/vector.cc
@@ -0,0 +1,20 @@
+#include "vector.hh"
+#include "string.hh"
+
+Vector::Vector(const Vector&n)
+          :dat(n.dat)
+{
+}    
+
+Vector
+Vector::operator-() const
+{
+    Vector v(*this); v*=-1; return v;
+}
+
+void
+Vector::set_unit(int j)
+{
+    fill(0.0);
+    dat[j] = 1.0;
+}
diff --git a/flower/vector.hh b/flower/vector.hh
new file mode 100644
index 0000000000..1929c674a9
--- /dev/null
+++ b/flower/vector.hh
@@ -0,0 +1,111 @@
+#ifndef VECTOR_HH
+#define VECTOR_HH
+
+#include <math.h>
+#include "real.hh"
+#include "vray.hh"
+
+class Dstream;
+class String;
+void set_matrix_debug(Dstream&ds);
+
+/// a row of numbers
+class Vector  {
+    svec<Real> dat;
+public:
+    void OK() const { dat.OK();}
+    int dim() const { return dat.sz(); }
+    Vector() { }
+    Vector(const Vector&n);
+    Vector(int n) {
+	dat.set_size(n);
+	fill(0);
+    }
+    void insert(Real v, int i) {
+	dat.insert(v,i);
+    }
+    void del(int i) { dat.del(i); }
+    operator String() const;
+    void fill(Real r) {
+	for (int i=0; i < dim(); i++)
+	    dat[i] =r;
+    }
+
+    void operator +=(Vector v) {
+	assert(v.dim() == dim());
+	for (int i=0; i < dim(); i++)
+	    dat[i] += v.dat[i];
+    }
+    
+    void operator /=(Real a) {
+	(*this) *= 1/a;
+    }
+
+    void operator *=(Real a) {
+	for (int i=0; i < dim(); i++)
+	    dat[i] *= a;
+    }
+
+    void operator -=(Vector v) {
+	assert(v.dim() == dim());
+	for (int i=0; i < dim(); i++)
+	    dat[i] -= v(i);	
+    }
+
+    Real &operator()(int i) { return dat[i]; }
+    Real operator()(int i) const { return dat[i]; }
+    Real elem(int i) { return dat[i]; }
+    Real operator *(Vector v) const {
+	Real ip=0;
+	assert(v.dim() == dim());
+	for (int i=0; i < dim(); i++)
+	    ip += dat[i] *v(i);
+	return ip;
+    }
+    Vector operator-() const;
+    Real norm() {
+	return sqrt(norm_sq() );
+    }
+    Real norm_sq() {
+	return ((*this) * (*this));
+    }
+    operator svec<Real> () { return dat; }
+    void print() const;
+    /// set to j-th element of unit-base
+    void set_unit(int j) ;
+};
+/**
+    a vector. Storage is handled in svec, Vector only does the mathematics.
+ */
+
+inline Vector
+operator+(Vector a, Vector const &b) {
+    a += b;
+    return a;
+}
+
+inline Vector
+operator-(Vector a, Vector const &b) {
+    a -= b;
+    return a;
+}
+
+inline Vector
+operator*(Vector v, Real a) {
+    v *= a;
+    return v;
+}
+
+inline Vector
+operator*( Real a,Vector v) {
+    v *= a;
+    return v;
+}
+
+inline Vector
+operator/(Vector v,Real a) {
+    v *= 1/a;
+    return v;
+}
+
+#endif
diff --git a/flower/vsmat.hh b/flower/vsmat.hh
new file mode 100644
index 0000000000..4e1ea797b6
--- /dev/null
+++ b/flower/vsmat.hh
@@ -0,0 +1,141 @@
+#ifndef VSMAT_HH
+#define VSMAT_HH
+#include "vray.hh"
+#include "real.hh"
+/// a  matrix storage baseclass.
+class virtual_smat {
+    
+
+public:
+    /// check invariants
+    virtual void OK() const=0;
+    
+    /// height of matrix
+    virtual int rows() const = 0;
+
+    /// width of matrix
+    virtual int cols() const = 0;
+    
+    /// set the size. contents lost
+    virtual void set_size(int i, int j) = 0;
+    /**      
+      PRE
+      i >=0, j>=0
+    */
+    
+    /// set the size to square dimen. contents lost
+    virtual void set_size(int i) = 0;
+    /**
+      PRE
+      i>=0
+    */
+    /// set the size to i
+    virtual void resize(int i, int j) = 0;
+    /**
+
+      keep contents. If enlarged contents unspecified
+        
+      PRE
+      i>=0, j>=0
+    
+    */
+
+    /// set the size to square dimen. contents kept
+    virtual void resize(int i) = 0;
+    /**    
+    Keep contents. If enlarged contents are unspecified
+    
+    PRE
+    i>=0  
+    */
+    
+    /// access an element
+    virtual Real& elem(int i,int j) = 0;
+    /**
+    access an element.
+
+    Generate an errormessage, if this happens
+    in the 0-part of a sparse matrix.
+    */
+
+    /// access a element, no modify
+    virtual const Real& elem(int i, int j) const = 0;
+
+#if 1
+    virtual svec<Real> row(int i) const = 0;
+    virtual svec<Real> column(int j) const = 0;
+#endif
+
+    /// add a row
+    virtual void insert_row(int k)=0;
+    /**
+    add a row to the matrix before  row k. Contents
+    of added row are unspecified
+
+      0 <= k <= rows()
+    */
+
+    /// delete a row 
+    virtual void delete_row(int k)=0;
+      /**
+      delete a row from this matrix.
+
+      PRE
+      0 <= k < rows();
+    */
+        virtual void delete_column(int k)=0;
+    virtual ~virtual_smat() { }
+    virtual virtual_smat *clone()=0;
+
+
+    /// is there a next?
+    virtual bool mult_ok(int i, int j) const=0;
+    /**
+      at end of matrix? when doing loop
+
+      for(i=0; i<h; i++)
+        for(j=0; j<w; j++)
+          ...
+
+    */
+    /// iterate
+    virtual void mult_next(int &i, int &j) const  = 0;
+    /**
+      walk through matrix (regular multiply)
+      get next j for row i, or get next row i and reset j.
+      this will make sparse matrix implementation easy.
+    
+      PRE
+      mult_ok(i,j)
+     */
+    virtual bool trans_ok(int i, int j) const=0;
+    /**
+      valid matrix entry. return false if at end of row
+    */
+    virtual void trans_next(int &i, int &j) const  = 0;
+    /**
+      walk through matrix (transposed multiply).
+      Get next i (for column j)
+    
+      PRE
+      ver_ok(i,j)
+     */
+
+    /// generate a "Full_storage" matrix    
+    static virtual_smat *get_full(int n, int m);
+
+};
+    
+/** base class for interface with matrix storageclasses.  There are no
+    iterators for matrixclasses, since matrices are (like arrays)
+    explicitly int-indexed.
+
+    Iteration is provided by *_next, *_ok, which update and check both
+    index variables simultaneously.
+
+    TODO
+    determine type of product matrix.
+
+*/
+
+#endif