*** empty log message ***

auxstorage.h

@@ -14,6 +14,8 @@
 //CAUTION:
 //it will not work if T is itself a class with dynamically allocated components
 //it cannot be implemented for SparseMat, which lacks fixed record length
+//for more complex I/O use put() and get() methods of the individual classes
+
 
 
 template<typename T>
@@ -38,6 +40,7 @@ public:
 template<typename T>
 AuxStorage<T>::AuxStorage(void)
 {
+//mkstemp probable does not support O_LARGEFILE?!
 strcpy(filename,"AUX_XXXXXX"); 
 mktemp(filename); 
 unlink(filename);

la_traits.h

@@ -1,40 +1,128 @@
 ////////////////////////////////////////////////////////////////////////////
-//traits classes
+//LA traits classes
 
 #ifndef _LA_TRAITS_INCL
 #define _LA_TRAITS_INCL
 
-//default one, good for numbers
+//forward declarations
-template<class C> struct NRMat_traits {
+template<typename C> class NRVec;
-typedef C elementtype;
+template<typename C> class NRMat;
-typedef C producttype;
+template<typename C> class NRSMat;
-static C norm (const  C &x) {return abs(x);}
+template<typename C> class SparseMat;
-static void axpy (C &s, const C &x, const C &c) {s+=x*c;}
+
+//let's do some simple template metaprogramming and preprocessing
+//to keep the thing general and compact
+
+typedef class scalar_false {};
+typedef class scalar_true {};
+
+//default is non-scalar
+template<typename C>
+class isscalar {
+	typedef scalar_false scalar_type;
 	};
 
 //specializations
-template<> struct NRMat_traits<NRMat<double> > {
+#define SCALAR(X) \
-typedef double elementtype;
+class isscalar<X> {typedef scalar_true scalar_type;};
-typedef NRMat<double> producttype;
+
-static double norm (const NRMat<double> &x) {return x.norm();}
+//declare what is scalar
-static void axpy (NRMat<double>&s, const NRMat<double> &x, const double c) {s.axpy(c,x);}
+SCALAR(char)
+SCALAR(short)
+SCALAR(int)
+SCALAR(long)
+SCALAR(long long)
+SCALAR(unsigned char)
+SCALAR(unsigned short)
+SCALAR(unsigned int)
+SCALAR(unsigned long)
+SCALAR(unsigned long long)
+SCALAR(float)
+SCALAR(double)
+SCALAR(complex<float>)
+SCALAR(complex<double>)
+SCALAR(void *)
+
+#undef SCALAR
+
+
+//now declare the traits for scalars and for composed classes
+template<typename C, typename Scalar> struct LA_traits_aux {};
+
+//complex scalars
+template<typename C>
+struct LA_traits_aux<complex<C>, scalar_true> {
+typedef complex<C> elementtype;
+typedef complex<C> producttype;
+typedef C normtype;
+static normtype norm (const  complex<C> &x) {return abs(x);}
+static void axpy (complex<C> &s, const complex<C> &x, const complex<C> &c) {s+=x*c;}
+static void get(int fd, complex<C> &x, bool dimensions=0) {if(sizeof(complex<C>)!=read(fd,&x,sizeof(complex<C>))) laerror("read error");}
+static void put(int fd, const complex<C> &x, bool dimensions=0) {if(sizeof(complex<C>)!=write(fd,&x,sizeof(complex<C>))) laerror("write error");}
+static void multiget(unsigned int n,int fd, complex<C> *x, bool dimensions=0){if((ssize_t)(n*sizeof(complex<C>))!=read(fd,x,n*sizeof(complex<C>))) laerror("read error");}
+static void multiput(unsigned int n, int fd, const complex<C> *x, bool dimensions=0) {if((ssize_t)(n*sizeof(complex<C>))!=write(fd,x,n*sizeof(complex<C>))) laerror("write error");}
+
 };
 
-template<> struct NRMat_traits<NRSMat<double> > {
+//non-complex scalars
-typedef double elementtype;
+template<typename C>
-typedef NRMat<double> producttype;
+struct LA_traits_aux<C, scalar_true> {
-static const double norm (const NRSMat<double> &x) {return x.norm(0.);}
+typedef C elementtype;
-static void axpy (NRSMat<double>&s, const NRSMat<double> &x, const double c) {s.axpy(c,x);}
+typedef C producttype;
+typedef C normtype;
+static normtype norm (const  C &x) {return abs(x);}
+static void axpy (C &s, const C &x, const C &c) {s+=x*c;}
+static void put(int fd, const C &x, bool dimensions=0) {if(sizeof(C)!=write(fd,&x,sizeof(C))) laerror("write error");}
+static void get(int fd, C &x, bool dimensions=0) {if(sizeof(C)!=read(fd,&x,sizeof(C))) laerror("read error");}
+static void multiput(unsigned int n,int fd, const C *x, bool dimensions=0){if((ssize_t)(n*sizeof(C))!=write(fd,x,n*sizeof(C))) laerror("write error");}
+static void multiget(unsigned int n, int fd, C *x, bool dimensions=0) {if((ssize_t)(n*sizeof(C))!=read(fd,x,n*sizeof(C))) laerror("read error");}
 };
 
 
-template<> struct NRMat_traits<NRMat<complex<double> > > {
+//prepare for non-scalar classes
-typedef complex<double> elementtype;
+
-typedef NRMat<complex<double> > producttype;
+template<typename C>
-static double norm (const NRMat<complex<double> >  &x) {return x.norm();}
+struct LA_traits; //forward declaration needed for template recursion
-static void axpy (NRMat<complex<double> >&s, const NRMat<complex<double> > &x, const complex<double> c) {s.axpy(c,x);}
+
+#define generate_traits(X) \
+template<typename C> \
+struct LA_traits_aux<X<C>, scalar_false> { \
+typedef C elementtype; \
+typedef X<C> producttype; \
+typedef typename LA_traits<C>::normtype normtype; \
+static normtype norm (const X<C> &x) {return x.norm();} \
+static 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) {x.put(fd,dimensions);} \
+static void get(int fd, C &x, bool dimensions=1) {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);} \
 };
 
 
+//non-scalar types defined in this library
+generate_traits(NRMat)
+generate_traits(NRVec)
+generate_traits(SparseMat)
+
+#undef generate_traits
+
+//non-scalar exceptions (smat product type)
+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 normtype norm (const NRSMat<C> &x) {return x.norm();} 
+static 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) {x.put(fd,dimensions);} 
+static void get(int fd, C &x, bool dimensions=1) {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);} \
+};
+
+
+//the final traits class
+template<typename C>
+struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {};
 
 #endif

laerror.cc

@@ -1,16 +1,30 @@
-// LA errorr handler
+// LA and general error handler
 #include <iostream>
 #include "laerror.h"
+#include <stdio.h>
+#include <errno.h>
 void laerror(const char *s1)
 {
   std::cerr << "LA:ERROR - ";
-  if(!s1)
+  std::cout << "LA:ERROR - ";
-    std::cerr << "udefined.\n";
+  if(s1)
-  else {
+  {
-    if(s1) std::cerr << s1;
+    std::cerr << s1  << "\n";
-	std::cerr << "\n";
+    std::cout << s1  << "\n";
   }
+  if(errno) perror("system error");
   throw LAerror(s1);
   exit(1);
 }
 
+//stub for f77 blas called from strassen routine
+extern "C" void xerbla_(const char name[6], int *n)
+{
+char msg[128];
+strcpy(msg,"LAPACK or BLAS error in routine ");
+strncat(msg,name,6);
+sprintf(msg+strlen(msg),": illegal value of parameter #%d",*n);
+laerror(msg);
+}
+
+

mat.cc

@@ -1,4 +1,12 @@
 #include "mat.h"
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+extern "C" {
+extern ssize_t read(int, void *, size_t);
+extern ssize_t write(int, const void *, size_t);
+}
 // TODO :
 //
 
@@ -14,6 +22,48 @@ template NRMat<char>;
  *  Templates first, specializations for BLAS next
  */
 
+//raw I/O
+template <typename T>
+void NRMat<T>::put(int fd, bool dim) const
+{
+errno=0;
+if(dim)
+{
+if(sizeof(int) != write(fd,&nn,sizeof(int))) laerror("cannot write");
+if(sizeof(int) != write(fd,&mm,sizeof(int))) laerror("cannot write");
+}
+LA_traits<T>::multiput(nn*mm,fd,
+#ifdef MATPTR
+        v[0]
+#else
+        v
+#endif
+,dim);
+}
+
+template <typename T>
+void NRMat<T>::get(int fd, bool dim)
+{
+int nn0,mm0;
+errno=0;
+if(dim)
+{
+if(sizeof(int) != read(fd,&nn0,sizeof(int))) laerror("cannot read");
+if(sizeof(int) != read(fd,&mm0,sizeof(int))) laerror("cannot read");
+resize(nn0,mm0);
+}
+else
+copyonwrite();
+LA_traits<T>::multiget(nn*mm,fd,
+#ifdef MATPTR
+        v[0]
+#else
+        v
+#endif
+,dim);
+}
+
+
 
 
 // Assign diagonal
@@ -817,7 +867,7 @@ istream& operator>>(istream  &s, NRMat<T> &x)
                 return s;
                 }
 
-
+//direct sum and product (oplus, otimes) to be done
 
 
 

mat.h

@@ -3,6 +3,7 @@
 
 #include "vec.h"
 #include "smat.h"
+#include "la_traits.h"
 
 template <typename T>
 class NRMat {
@@ -29,6 +30,12 @@ public:
 	NRMat(const NRVec<T> &rhs, const int n, const int m);
 #endif
 	~NRMat();
+#ifdef MATPTR
+	const bool operator!=(const NRMat &rhs) const {if(nn!=rhs.nn || mm!=rhs.mm) return 1; return(memcmp(v[0],rhs.v[0],nn*mm*sizeof(T)));}
+#else
+        const bool operator!=(const NRMat &rhs) const {if(nn!=rhs.nn || mm!=rhs.mm) return 1; return(memcmp(v,rhs.v,nn*mm*sizeof(T)));}
+#endif
+	const bool operator==(const NRMat &rhs) const {return !(*this != rhs);};
 	inline int getcount() const {return count?*count:0;}
 	NRMat & operator=(const NRMat &rhs);  //assignment
 	NRMat & operator=(const T &a);    //assign a to diagonal
@@ -50,6 +57,8 @@ public:
 	inline const NRMat operator-(const NRSMat<T> &rhs) const;
 	const T dot(const NRMat &rhs) const; // scalar product of Mat.Mat
 	const NRMat operator*(const NRMat &rhs) const; // Mat * Mat
+	const NRMat oplus(const NRMat &rhs) const; //direct sum
+	const NRMat otimes(const NRMat &rhs) const; //direct product
 	void diagmultl(const NRVec<T> &rhs); //multiply by a diagonal matrix from L
 	void diagmultr(const NRVec<T> &rhs); //multiply by a diagonal matrix from R
 	const NRMat operator*(const NRSMat<T> &rhs) const; // Mat * Smat
@@ -65,6 +74,8 @@ public:
 	inline const T& operator()(const int i, const int j) const;
 	inline int nrows() const;
 	inline int ncols() const;
+        void get(int fd, bool dimensions=1);
+        void put(int fd, bool dimensions=1) const;
 	void copyonwrite();
 	void resize(const int n, const int m);
 	inline operator T*(); //get a pointer to the data
@@ -446,9 +457,24 @@ template <typename T>
 void NRMat<T>::resize(const int n, const int m)
 {
 #ifdef DEBUG
-        if (n<=0 || m<=0) laerror("illegal dimensions in Mat::resize()");
+        if (n<0 || m<0 || n>0 && m==0 || n==0 && m>0) laerror("illegal dimensions in Mat::resize()");
 #endif
         if (count)
+	    {
+	        if(n==0 && m==0)
+        		{
+		        if(--(*count) <= 0) {
+#ifdef MATPTR
+                		if(v) delete[] (v[0]);
+#endif
+                		if(v) delete[] v;
+		                delete count;
+		                }
+		        count=0;
+		        nn=mm=0;
+		        v=0;
+		        return;
+		        }
                 if (*count > 1) {
                         (*count)--;
                         count = 0;
@@ -456,6 +482,7 @@ void NRMat<T>::resize(const int n, const int m)
                         nn = 0;
                         mm = 0;
                 }
+	    }
         if (!count) {
                 count = new int;
                 *count = 1;

matexp.h

@@ -6,7 +6,6 @@
 // is defined containing definition of an element type, norm and axpy operation
 
 #include "la_traits.h"
-#include "sparsemat_traits.h"
 
 template<class T,class R>
 const T polynom2(const T &x, const NRVec<R> &c)
@@ -68,7 +67,7 @@ for(i=0; i<=n/m;i++)
 		if(k>n) break;
 		if(j==0) {if(i==0) s=x; /*just to get the dimensions of the matrix*/ s=c[k]; /*create diagonal matrix*/}
 		else  
-			NRMat_traits<T>::axpy(s,xpows[j-1],c[k]); //general  s+=xpows[j-1]*c[k]; but more efficient for matrices
+			LA_traits<T>::axpy(s,xpows[j-1],c[k]); //general  s+=xpows[j-1]*c[k]; but more efficient for matrices
 		}
 
 	if(i==0) {r=s; f=xpows[m-1];}
@@ -125,7 +124,7 @@ template<class T>
 const T ipow( const T &x, int i)
 {
 if(i<0) laerror("negative exponent in ipow");
-if(i==0) {T r=x; r=(T)1; return r;}//trick for matrix dimension
+if(i==0) {T r=x; r=(typename LA_traits<T>::elementtype)1; return r;}//trick for matrix dimension
 if(i==1) return x;
 T y,z;
 z=x;
@@ -153,7 +152,7 @@ return int(ceil(log(n)/log2-log(.75)));
 
 
 template<class T>
-NRVec<typename NRMat_traits<T>::elementtype> exp_aux(const T &x, int &power)
+NRVec<typename LA_traits<T>::elementtype> exp_aux(const T &x, int &power)
 {
 //should better be computed by mathematica to have accurate last digits, chebyshev instead, see exp in glibc
 static double exptaylor[]={
@@ -179,7 +178,7 @@ static double exptaylor[]={
 8.2206352466243294955e-18,
 4.1103176233121648441e-19,
 0.};
-double mnorm= NRMat_traits<T>::norm(x);
+double mnorm= LA_traits<T>::norm(x);
 power=nextpow2(mnorm);
 double scale=exp(-log(2.)*power);
 
@@ -198,7 +197,7 @@ while(t*exptaylor[n]>precision);//taylor 0 will terminate in any case
 
 
 int i; //adjust the coefficients in order to avoid scaling the argument
-NRVec<typename NRMat_traits<T>::elementtype> taylor2(n+1);
+NRVec<typename LA_traits<T>::elementtype> taylor2(n+1);
 for(i=0,t=1.;i<=n;i++)
 	{
 	taylor2[i]=exptaylor[i]*t;
@@ -215,7 +214,7 @@ const T exp(const T &x)
 int power;
 
 //prepare the polynom of and effectively scale T
-NRVec<typename NRMat_traits<T>::elementtype> taylor2=exp_aux(x,power);
+NRVec<typename LA_traits<T>::elementtype> taylor2=exp_aux(x,power);
 
 T r=polynom(x,taylor2); //for accuracy summing from the smallest terms up would be better, but this is more efficient for matrices
 
@@ -233,7 +232,7 @@ const V exptimes(const M &mat, V vec) //uses just matrix vector multiplication
 if(mat.nrows()!=mat.ncols()||(unsigned int) mat.nrows() != (unsigned int)vec.size()) laerror("inappropriate sizes in exptimes");
 int power;
 //prepare the polynom of and effectively scale the matrix
-NRVec<typename NRMat_traits<M>::elementtype> taylor2=exp_aux(mat,power);
+NRVec<typename LA_traits<M>::elementtype> taylor2=exp_aux(mat,power);
 
 V result(mat.nrows());
 for(int i=1; i<=(1<<power); ++i) //unfortunatelly, here we have to repeat it many times, unlike if the matrix is stored explicitly

nonclass.h

@@ -107,9 +107,9 @@ const NRMat<T> inverse(NRMat<T> a, T *det=0)
 
 //general determinant
 template<class MAT>
-const typename NRMat_traits<MAT>::elementtype determinant(MAT a)//again passed by value
+const typename LA_traits<MAT>::elementtype determinant(MAT a)//again passed by value
 {
-typename NRMat_traits<MAT>::elementtype det;
+typename LA_traits<MAT>::elementtype det;
 if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
 linear_solve(a,NULL,&det);
 return det;

smat.cc

@@ -1,4 +1,12 @@
 #include "smat.h"
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+extern "C" {
+extern ssize_t read(int, void *, size_t);
+extern ssize_t write(int, const void *, size_t);
+}
 // TODO
 // specialize unary minus
 
@@ -17,6 +25,35 @@ template NRSMat<char>;
  *
  */
 
+//raw I/O
+template <typename T>
+void NRSMat<T>::put(int fd, bool dim) const
+{
+errno=0;
+if(dim)
+{
+if(sizeof(int) != write(fd,&nn,sizeof(int))) laerror("cannot write");
+if(sizeof(int) != write(fd,&nn,sizeof(int))) laerror("cannot write");
+}
+LA_traits<T>::multiput(NN2,fd,v,dim);
+}
+
+template <typename T>
+void NRSMat<T>::get(int fd, bool dim)
+{
+int nn0[2]; //align at least 8-byte
+errno=0;
+if(dim)
+{
+if(2*sizeof(int) != read(fd,&nn0,2*sizeof(int))) laerror("cannot read");
+resize(nn0[0]);
+}
+else 
+copyonwrite();
+LA_traits<T>::multiget(NN2,fd,v,dim);
+}
+
+
 // conversion ctor, symmetrize general Mat into SMat
 template <typename T>
 NRSMat<T>::NRSMat(const NRMat<T> &rhs)

smat.h

@@ -3,6 +3,7 @@
 
 #include "vec.h"
 #include "mat.h"
+#include "la_traits.h"
 
 #define NN2 (nn*(nn+1)/2)
 template <class T>
@@ -25,6 +26,8 @@ public:
 	NRSMat & operator|=(const NRSMat &rhs);	//assignment to a new copy
 	NRSMat & operator=(const NRSMat &rhs);	//assignment
 	NRSMat & operator=(const T &a);		//assign a to diagonal
+        const bool operator!=(const NRSMat &rhs) const {if(nn!=rhs.nn) return 1; return(memcmp(v,rhs.v,NN2*sizeof(T)));}
+        const bool operator==(const NRSMat &rhs) const {return !(*this != rhs);};
 	inline NRSMat & operator*=(const T &a);
 	inline NRSMat & operator+=(const T &a); 
 	inline NRSMat & operator-=(const T &a); 
@@ -54,6 +57,8 @@ public:
 	void axpy(const T alpha, const NRSMat &x); // this+= a*x
 	inline const T amax() const;
 	const T trace() const;
+	void get(int fd, bool dimensions=1);
+        void put(int fd, bool dimensions=1) const;
 	void copyonwrite();
 	void resize(const int n);
 	inline operator T*(); //get a pointer to the data
@@ -396,15 +401,29 @@ template <typename T>
 void NRSMat<T>::resize(const int n)
 {
 #ifdef DEBUG
-        if (n <= 0) laerror("illegal matrix dimension in resize of Smat");
+        if (n < 0) laerror("illegal matrix dimension in resize of Smat");
 #endif
         if (count)
+	{
+	    	if(n==0)
+	        {
+	        if(--(*count) <= 0) {
+	                if(v) delete[] (v);
+	                delete count;
+	                }
+	        count=0;
+	        nn=0;
+	        v=0;
+	        return;
+	        }
+
                 if(*count > 1) {        //detach from previous
                         (*count)--;
                         count = 0;
                         v = 0;
                         nn = 0;
                 }
+	}
         if (!count)     {                               //new uninitialized vector or just detached
                 count = new int;
                 *count = 1;

sparsemat.cc

@@ -2,7 +2,10 @@
 #include <cmath>
 #include <complex>
 #include <iostream>
-
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
 #include "sparsemat.h"
 
 //////////////////////////////////////////////////////////////////////////////
@@ -56,6 +59,72 @@ istream& operator>>(istream  &s, SparseMat<T> &x)
                 return s;
                 }
 
+extern "C" {
+extern ssize_t read(int, void *, size_t);
+extern ssize_t write(int, const void *, size_t);
+}
+
+
+export template <class T>
+void SparseMat<T>::get(int fd, bool dimen)
+{
+errno=0;
+SPMatindex dim[2];
+if(dimen)
+{
+if(2*sizeof(SPMatindex) != read(fd,&dim,2*sizeof(SPMatindex))) laerror("cannot read");
+resize(dim[0],dim[1]);
+int symnon[2];
+if(2*sizeof(int) != read(fd,&symnon,2*sizeof(int))) laerror("cannot read");
+symmetric=symnon[0];
+nonzero=symnon[1];
+}
+else
+copyonwrite();
+
+matel<T> *l=NULL;
+do
+          {
+	  if(2*sizeof(SPMatindex) != read(fd,&dim,2*sizeof(SPMatindex))) laerror("cannot read");
+	  if(dim[0]+1==0 && dim[1]+1==0) break;
+          matel<T> *ll = l;
+          l= new matel<T>;
+          l->next= ll;
+          l->row=dim[0];
+          l->col=dim[1];
+	  LA_traits<T>::get(fd,l->elem,dimen); //general way to work when elem is some complex class again
+          } while(1);
+list=l;
+}
+
+
+export template <class T>
+void SparseMat<T>::put(int fd,bool dimen) const
+{
+errno=0;
+if(dimen)
+{
+if(sizeof(SPMatindex) != write(fd,&nn,sizeof(SPMatindex))) laerror("cannot write");
+if(sizeof(SPMatindex) != write(fd,&mm,sizeof(SPMatindex))) laerror("cannot write");
+int symnon[2];
+symnon[0]=symmetric;
+symnon[1]=nonzero;
+if(2*sizeof(int) != write(fd,symnon,2*sizeof(int))) laerror("cannot write");
+}
+matel<T> *l=list;
+while(l)
+	{
+	if(sizeof(SPMatindex) != write(fd,&l->row,sizeof(SPMatindex))) laerror("cannot write");
+	if(sizeof(SPMatindex) != write(fd,&l->col,sizeof(SPMatindex))) laerror("cannot write");
+	LA_traits<T>::put(fd,l->elem,dimen);//general way to work when elem is some non-scalar class again
+	l=l->next;
+	}
+SPMatindex sentinel[2];
+sentinel[0]=sentinel[1]=(SPMatindex) -1;
+if(2*sizeof(SPMatindex) != write(fd,&sentinel,2*sizeof(SPMatindex))) laerror("cannot write");
+}
+
+
 //helpers to be used from different functions
 export template <class T>
 void SparseMat<T>::unsort()
@@ -314,21 +383,35 @@ unsort(); //since there were NULLs introduced, rowsorted is not dense
 export template <class T>
 void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m)
 {
-        if(n<=0 || m<=0) laerror("illegal matrix dimension");
 	unsort(); 
         if(count) 
 		{
 		if(*count > 1) {(*count)--; count=NULL; list=NULL;} //detach from previous
 		else if(*count==1) deletelist();
+		if(count) delete count;
 		}
         nn=n;
         mm=m;
-        count=new int(1); //empty but defined matrix
+        if(nn||mm) count=new int(1); //empty but defined matrix
 	list=NULL;
 	symmetric=0;
+	nonzero=0;
 	colsorted=rowsorted=NULL;
 }
 
+export template <class T>
+void SparseMat<T>::incsize(const SPMatindex n, const SPMatindex m)
+{
+	if(symmetric && n!=m) laerror("unsymmetric size increment of a symmetric sparsemat");
+	if(!count && nn==0 && mm==0) count=new int(1);
+	copyonwrite();//this errors if !count 
+	unsort();
+        nn+=n;
+        mm+=m;
+}
+
+
+
 export template <class T>
 void SparseMat<T>::addsafe(const SPMatindex n, const SPMatindex m, const T elem)
 {
@@ -1072,17 +1155,80 @@ for(i=0; i<na;i++)
 simplify();
 }
 
+//direct sum and product -- only partly implemented at the moment
+export template<typename T>
+SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs)
+{
+if(issymmetric()) laerror("oplusequal symmetric-unsymmetric");
+SPMatindex n0=nn;
+SPMatindex m0=mm;
+incsize(rhs.nrows(), rhs.ncols());
+T tmp;
+for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
+	for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
+#ifdef SPARSEEPSILON
+		if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
+#else
+		if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
+#endif
+return *this;
+}
 
+export template<typename T>
+SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs)
+{
+if(!issymmetric()) laerror("oplusequal symmetric-unsymmetric");
+SPMatindex n0=nn;
+SPMatindex m0=mm;
+T tmp;
+incsize(rhs.nrows(), rhs.ncols());
+for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
+        for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
+#ifdef SPARSEEPSILON
+                if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
+#else
+                if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
+#endif
+return *this;
+}
+
+export template <class T>
+SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs)
+{
+if(symmetric != rhs.symmetric) laerror("incompatible symmetry of sparsemats in oplusequal");
+if(!count) {count=new int;  *count=1; list=NULL;}
+
+SPMatindex n0=nn;
+SPMatindex m0=mm; 
+incsize(rhs.nrows(), rhs.ncols());
+
+register matel<T> *l=rhs.list;
+while(l)
+	{
+#ifdef SPARSEEPSILON
+	if(abs(l->elem)>SPARSEEPSILON) 
+#endif
+		add(n0+l->row,m0+l->col,l->elem);
+	l=l->next;
+	}
+return *this;
+}
 
 
 #ifdef _GLIBCPP_NO_TEMPLATE_EXPORT
 
 #define INSTANTIZE(T) \
+template SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs);\
+template SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs);\
+template SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs);\
 template ostream& operator<<(ostream &s, const SparseMat<T> &x); \
 template istream& operator>>(istream &s, SparseMat<T> &x); \
+template void SparseMat<T>::get(int fd, bool dimen); \
+template void SparseMat<T>::put(int fd, bool dimen) const; \
 template void SparseMat<T>::copyonwrite(); \
-template void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m); \
 template void SparseMat<T>::unsort(); \
+template void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m); \
+template void SparseMat<T>::incsize(const SPMatindex n, const SPMatindex m); \
 template unsigned int SparseMat<T>::sort(int type) const; \
 template unsigned int SparseMat<T>::length() const; \
 template void SparseMat<T>::deletelist(); \

sparsemat.h

@@ -90,6 +90,15 @@ public:
 	inline const NRVec<T> operator*(const NRVec<T> &rhs) const {return multiplyvector(rhs);} //sparse matrix * dense vector
 	void diagonalof(NRVec<T> &, const bool divide=0) const; //get diagonal
 	const SparseMat operator*(const SparseMat &rhs) const; 
+	SparseMat & oplusequal(const SparseMat &rhs); //direct sum
+	SparseMat & oplusequal(const NRMat<T> &rhs);
+	SparseMat & oplusequal(const NRSMat<T> &rhs);
+	const SparseMat oplus(const SparseMat &rhs) const {return SparseMat(*this).oplusequal(rhs);}; //direct sum
+	const SparseMat oplus(const NRMat<T> &rhs) const {return SparseMat(*this).oplusequal(rhs);};
+	const SparseMat oplus(const NRSMat<T> &rhs) const {return SparseMat(*this).oplusequal(rhs);};
+	const SparseMat otimes(const SparseMat &rhs) const; //direct product
+        const SparseMat otimes(const NRMat<T> &rhs) const;
+        const SparseMat otimes(const NRSMat<T> &rhs) const;
         void gemm(const T beta, const SparseMat &a, const char transa, const SparseMat &b, const char transb, const T alpha);//this := alpha*op( A )*op( B ) + beta*this, if this is symemtric, only half will be added onto it
 	const T dot(const SparseMat &rhs) const; //supervector dot product
 	const T dot(const NRMat<T> &rhs) const; //supervector dot product
@@ -99,7 +108,10 @@ public:
 	void setlist(matel<T> *l) {list=l;}
 	inline SPMatindex nrows() const {return nn;}
         inline SPMatindex ncols() const {return mm;}
-	void resize(const SPMatindex n, const SPMatindex m);
+	void get(int fd, bool dimensions=1);
+	void put(int fd, bool dimensions=1) const;
+	void resize(const SPMatindex n, const SPMatindex m); //destructive
+	void incsize(const SPMatindex n, const SPMatindex m); //increase size without destroying the data
 	void transposeme();
 	const SparseMat transpose() const;
 	inline void setsymmetric() {if(nn!=mm) laerror("non-square cannot be symmetric"); symmetric=1;}

strassen.cc

@@ -22,10 +22,3 @@ copyonwrite();
 fmm(transb,transa,mm,nn,k,alpha,b,b.mm, a, a.mm, beta,*this, mm,NULL,0);
 }
 
-//stub for f77 blas called from strassen routine
-extern "C" void xerbla_(const char *msg)
-{
-laerror(msg);
-}
-
-

vec.cc

@@ -1,5 +1,14 @@
 #include <iostream>
 #include "vec.h"
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+extern "C" {
+extern ssize_t read(int, void *, size_t);
+extern ssize_t write(int, const void *, size_t);
+}
+
 
 //////////////////////////////////////////////////////////////////////////////
 //// forced instantization in the corespoding object file
@@ -10,7 +19,11 @@ template istream & operator>>(istream  &s, NRVec< T > &x); \
 INSTANTIZE(double)
 INSTANTIZE(complex<double>)
 INSTANTIZE(int)
+INSTANTIZE(unsigned int)
+INSTANTIZE(short)
+INSTANTIZE(unsigned short)
 INSTANTIZE(char)
+INSTANTIZE(unsigned char)
 template NRVec<double>;
 template NRVec<complex<double> >;
 template NRVec<int>;
@@ -36,7 +49,7 @@ NRVec<T>::NRVec(const NRMat<T> &rhs)
 
 
 
-// ostream << NRVec
+// formatted I/O
 template <typename T>
 ostream & operator<<(ostream &s, const NRVec<T> &x)
 {
@@ -48,7 +61,6 @@ ostream & operator<<(ostream &s, const NRVec<T> &x)
   return s;
 }
 
-// istream >> NRVec
 template <typename T>
 istream & operator>>(istream &s, NRVec<T> &x)
 {
@@ -60,6 +72,39 @@ istream & operator>>(istream &s, NRVec<T> &x)
   return s;
 }
 
+
+//raw I/O
+template <typename T>
+void NRVec<T>::put(int fd, bool dim) const
+{
+errno=0;
+int pad=1; //align at least 8-byte
+if(dim)
+{
+if(sizeof(int) != write(fd,&nn,sizeof(int))) laerror("cannot write");
+if(sizeof(int) != write(fd,&pad,sizeof(int))) laerror("cannot write");
+}
+LA_traits<T>::multiput(nn,fd,v,dim);
+}
+
+template <typename T>
+void NRVec<T>::get(int fd, bool dim)
+{
+int nn0[2]; //align at least 8-byte
+errno=0;
+if(dim)
+{
+if(2*sizeof(int) != read(fd,&nn0,2*sizeof(int))) laerror("cannot read");
+resize(nn0[0]);
+}
+else
+copyonwrite();
+LA_traits<T>::multiget(nn,fd,v,dim);
+}
+
+
+
+
 // formatted print for NRVec
 template<typename T>
 void NRVec<T>::fprintf(FILE *file, const char *format, const int modulo) const
@@ -68,7 +113,7 @@ void NRVec<T>::fprintf(FILE *file, const char *format, const int modulo) const
 }
 
 // formatted scan for NRVec
-template <class T>
+template <typename T>
 void NRVec<T>::fscanf(FILE *f, const char *format)
 {
 	int n;

vec.h

@@ -2,7 +2,6 @@
 #define _LA_VEC_H_
 
 #include "laerror.h"
-
 extern "C" {
 #include "cblas.h"
 }
@@ -13,6 +12,8 @@ extern "C" {
 
 using namespace std;
 
+#include "la_traits.h"
+
 template <typename T> class NRVec;
 template <typename T> class NRSMat;
 template <typename T> class NRMat;
@@ -69,6 +70,8 @@ public:
 	NRVec & operator=(const NRVec &rhs);
 	NRVec & operator=(const T &a);  //assign a to every element
 	NRVec & operator|=(const NRVec &rhs);
+	const bool operator!=(const NRVec &rhs) const {if(nn!=rhs.nn) return 1; return(memcmp(v,rhs.v,nn*sizeof(T)));}
+	const bool operator==(const NRVec &rhs) const {return !(*this != rhs);};
 	const NRVec operator-() const;
 	inline NRVec & operator+=(const NRVec &rhs);
 	inline NRVec & operator-=(const NRVec &rhs);
@@ -99,6 +102,8 @@ public:
 			const T alpha, const NRVec &x);
 	void copyonwrite();
 	void resize(const int n);
+        void get(int fd, bool dimensions=1);
+        void put(int fd, bool dimensions=1) const;
 	NRVec & normalize();
 	inline const double norm() const;
 	inline const T amax() const;
@@ -495,15 +500,28 @@ template <typename T>
 void NRVec<T>::resize(const int n)
 {
 #ifdef DEBUG
-  if(n<=0) laerror("illegal vector dimension");
+  if(n<0) laerror("illegal vector dimension");
 #endif
   if(count)
+  {
+    if(n==0)
+	{
+	if(--(*count) <= 0) {
+                if(v) delete[] (v);
+                delete count;
+        	}
+	count=0;
+	nn=0;
+	v=0;
+	return;
+	}
     if(*count > 1) {
       (*count)--;
       count = 0;
       v = 0;
       nn = 0;
     }
+  }
   if(!count) {
     count = new int;
     *count = 1;