diff --git a/diis.h b/diis.h
index d6f09b8..90b9666 100644
--- a/diis.h
+++ b/diis.h
@@ -9,7 +9,9 @@
 #include "la_traits.h"
 #include "auxstorage.h"
 
-// T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields
+// Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields
+// actually it can be anything what has operator=(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype
+// and get() and put() if external storage is requested
 
 template<typename T>
 class DIIS
@@ -23,7 +25,9 @@ class DIIS
 	AuxStorage<Te> *st;
 	T *stor;
 public:
+	DIIS() {dim=0;}; //for array of diis
 	DIIS(const int n, const bool core=1);
+	void setup(const int n, const bool core=1);
 	~DIIS();
 	typename LA_traits<T>::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm
 	};
@@ -35,9 +39,21 @@ st=incore?NULL: new AuxStorage<Te>;
 stor= incore? new T[dim] : NULL;
 bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1; 
 aktdim=cyclicshift=0;
-
 }
 
+template<typename T>
+void DIIS<T>::setup(const int n, const bool core)
+{
+dim=n;
+incore=core;
+bmat.resize(n);
+st=incore?NULL: new AuxStorage<Te>;
+stor= incore? new T[dim] : NULL;
+bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1;
+aktdim=cyclicshift=0;
+}
+
+
 template<typename T>
 DIIS<T>::~DIIS()
 {
@@ -49,6 +65,7 @@ if(stor) delete[] stor;
 template<typename T>
 typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec)
 {
+if(!dim) laerror("attempt to extrapolate from uninitialized DIIS");
 //if dim exceeded, shift 
 if(aktdim==dim)
 	{
@@ -68,12 +85,12 @@ if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000;
 
 //calculate difference; 
 vec.copyonwrite();
-if(incore) vec -= stor[(aktdim-2+cyclicshift)%dim];
+if(incore) vec.axpy((Te)-1,stor[(aktdim-2+cyclicshift)%dim]);
 else
 	{
 	T tmp=vec;
 	st->get(tmp,(aktdim-2+cyclicshift)%dim);
-	vec -= tmp;
+	vec.axpy((Te)-1,tmp);
 	}
 
 //calculate overlaps of differences (if storage is cheap, they could rather be stored than recomputed)
diff --git a/fourindex.h b/fourindex.h
index 4008494..fb11e75 100644
--- a/fourindex.h
+++ b/fourindex.h
@@ -600,8 +600,17 @@ istream& operator>>(istream  &s, fourindex<I,T> &x)
 //note - loops for the twoelectronrealmullikan integral to be unique and in canonical order
 // i=1..n, j=1..i, k=1..i, l=1..(i==k?j:k)
 
+//general template declaration
 template<fourindexsymtype S, class T, class DUMMY> class fourindex_dense;
 
+//traits class
+template<fourindexsymtype S, class T, class DUMMY>
+struct LA_traits<fourindex_dense<S,T,DUMMY> > {
+typedef T elementtype;
+typedef typename LA_traits<T>::normtype normtype;
+};
+
+
 //make it as a derived class in order to be able to use it in a base class context - "supermatrix" operations
 template<class T, class I> 
 class fourindex_dense<twoelectronrealmullikan,T,I> : public NRSMat<T> {
@@ -768,7 +777,4 @@ if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
 
 
 
-
-
-
 #endif /*_fourindex_included*/
diff --git a/la_traits.h b/la_traits.h
index 24651ba..e415f63 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -25,7 +25,9 @@ extern "C" {
 //forward declarations
 template<typename C> class NRVec;
 template<typename C> class NRMat;
+template<typename C> class NRMat_from1;
 template<typename C> class NRSMat;
+template<typename C> class NRSMat_from1;
 template<typename C> class SparseMat;
 
 //for general sortable classes
@@ -174,30 +176,34 @@ static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest
 
 //non-scalar types defined in this library
 generate_traits(NRMat)
+generate_traits(NRMat_from1)
 generate_traits(NRVec)
 generate_traits(SparseMat)
 
 #undef generate_traits
 
 //smat
-template<typename C> 
-struct LA_traits_aux<NRSMat<C>, scalar_false> { 
-typedef C elementtype; 
-typedef NRMat<C> producttype; 
-typedef typename LA_traits<C>::normtype normtype; 
-static bool gencmp(const C *x, const C *y, int n) {for(int i=0; i<n; ++i) if(x[i]!=y[i]) return true; return false;}
-static inline bool bigger(const  C &x, const C &y) {return x>y;}
-static inline bool smaller(const  C &x, const C &y) {return x<y;}
-static inline normtype norm (const NRSMat<C> &x) {return x.norm();} 
-static inline void axpy (NRSMat<C>&s, const NRSMat<C> &x, const C c) {s.axpy(c,x);} 
-static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);} 
-static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);} 
-static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].put(fd,dimensions);} 
-static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].get(fd,dimensions);} 
-static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} 
-static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();}
+#define generate_traits_smat(X) \
+template<typename C>  \
+struct LA_traits_aux<X<C>, scalar_false> {  \
+typedef C elementtype;  \
+typedef NRMat<C> producttype;  \
+typedef typename LA_traits<C>::normtype normtype;  \
+static bool gencmp(const C *x, const C *y, int n) {for(int i=0; i<n; ++i) if(x[i]!=y[i]) return true; return false;} \
+static inline bool bigger(const  C &x, const C &y) {return x>y;} \
+static inline bool smaller(const  C &x, const C &y) {return x<y;} \
+static inline normtype norm (const X<C> &x) {return x.norm();}  \
+static inline void axpy (X<C>&s, const X<C> &x, const C c) {s.axpy(c,x);}  \
+static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);}  \
+static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);}  \
+static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].put(fd,dimensions);}  \
+static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].get(fd,dimensions);}  \
+static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];}  \
+static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();} \
 };
 
+generate_traits_smat(NRSMat)
+generate_traits_smat(NRSMat_from1)
 
 //the final traits class
 template<typename C>