csrmat claude-coded

2026-11-05 23:23:55 +01:00
1 parent 90bd68b5ef
commit 74479c1a27
@@ -66,6 +66,12 @@ decisions and plans agreed with the author, so that they survive a fresh clone.
  adds x at (i,j) and conj(x) at (j,i); NRSMat(SparseSMat) takes the lower triangle; operator* is the
  general product A B (correct for hermitian, complex symmetric exp(iH) and general matrices);
  SparseSMat::gemm computes A^T B (A^\dagger B with 'c'); otimes_sparse = M + M^\dagger incl. doubled diagonal.
+- CSRMat (completed 2026-11): compressed row format, members NRVec rowptr/col/v (share independently, value-only operations keep
+  the index arrays shared), columns strictly ascending in each row, full matrix always stored; flag symmetric only DECLARES hermitian
+  (A=A^dagger, for davidson/lanczos), set by CSRMat(NRSMat) and for real types from a symmetric SparseMat, never checked. Built by
+  conversion (SparseMat, SparseSMat, NRMat, NRSMat, arrays); conversions drop |x|<=sparse_epsilon, arithmetic keeps cancelled zeros
+  (simplify). dot conjugates the left operand (as NRMat); permuteme_rows permutes only rows; CSR*CSR = Gustavson; types double, float,
+  complex<double>, complex<float>, CPU only (cuSPARSE format for later). Tests bugtests/core/csrmat1,3,4,5 and iterative/csrmat2.
 - SparseMat with the symmetric flag is complex symmetric (no conjugation); its permuteme_* members use
  0-based matrix indices, NRPerm/CyclePerm are always 1-based.
 - Permutations: applypermutation(p,v,inverse) uses the same convention as NRVec::permuted (r_{p(i)}=v_i
@@ -83,8 +89,7 @@ decisions and plans agreed with the author, so that they survive a fresh clone.
 - a-X and a/X (scalar minus/divided by object) are deleted on purpose; a single integer constructor is
  forbidden for NRMat, NRMat_from1, Mat3 (Mat3(double) works).
 - Real-storage gdiagonalize uses the dgeev convention (Re + i Im).
- The NRMat offset view constructor is a feature. MATPTR is kept (rarely used). CSRMat is unfinished and
-  documented as unsupported, but kept. NRMat::strassen (2026-11-06): native Strassen-Winograd in strassen.cc, gemm-compatible
+- The NRMat offset view constructor is a feature. MATPTR is kept (rarely used). CSRMat: see above. NRMat::strassen (2026-11-06): native Strassen-Winograd in strassen.cc, gemm-compatible
  parameters + cutoff (default LA_traits<T>::strassen_cutoff()=2048 all types, to be tuned on HPC), explicit only (gemm stays classical),
  CBLAS leaves on strided quadrants, device primitives separated for CUDA; NO_STRASSEN removed everywhere. CUDA (CUDALA) is deferred/experimental
  until the author has NVIDIA RTX hardware (known CUDA issues collected in suspected_bugs.md F10).
@@ -98,7 +103,7 @@ decisions and plans agreed with the author, so that they survive a fresh clone.
  conjgrad, gmres, davidson, lanczos, diis: no changes needed beyond the defaults). Float twins are generated from the double specializations (placed right after them). Precision-
  dependent defaults live in LA_precision<T> (la_traits.h, only double and float; double values unchanged). optimize done
  2026-11-05: own trait opt_precision<T> in optimize.h, INSTANTIZE(float), float tests in testoptimize.cc (test_float);
-  float LM on badly scaled problems needs my_scaling=true. reg/regsurf, csrmat, CUDA stay double-only.
+  float LM on badly scaled problems needs my_scaling=true. reg/regsurf, CUDA stay double-only.
  Tensor Tucker bugs fixed 2026-11-06 (complex: factor = vt.transpose(false); unnamed tensors with inverseorder; bugtests/tensor/tucker.cc).
  BLAS-based NRVec members are meant for floating types only, not integers.
 - NRSMat: possible bool member `hermitian` toggling hermitian vs complex symmetric behaviour
@@ -1,6 +1,6 @@
 /*
    LA: linear algebra C++ interface library
-    Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+    Copyright (C) 2008-2026 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -24,38 +24,955 @@
 #include <fcntl.h>
 #include <errno.h>
 #include "csrmat.h"
+#include "nonclass.h"
+#include <algorithm>

 namespace LA {

+//raw pointers of possibly empty arrays (NULL when empty)
+template<typename U> static inline const U *s_ptr(const NRVec<U> &x) {return x.size()? (const U *)x : NULL;}
+template<typename U> static inline U *s_wptr(NRVec<U> &x) {return x.size()? (U *)x : NULL;}
+
+template <typename T>
+CSRMat<T>::CSRMat(const SPMatindex n, const SPMatindex m)
+: nn(n), mm(m), symmetric(false), rowptr(n+1)
+{
+rowptr.clear();
+}
+
+template <typename T>
+CSRMat<T>::CSRMat(const SPMatindex n)
+: nn(n), mm(n), symmetric(false), rowptr(n+1)
+{
+rowptr.clear();
+}
+
+template <typename T>
+CSRMat<T>::CSRMat(const SPMatindex n, const SPMatindex m, const NRVec<SPMatindex> &rp, const NRVec<SPMatindex> &c, const NRVec<T> &x)
+: nn(n), mm(m), symmetric(false), rowptr(rp), col(c), v(x)
+{
+#ifdef DEBUG
+if(rp.size()!=(int)n+1 || rp[0]!=0 || (int)rp[n]!=c.size() || c.size()!=x.size()) laerror("inconsistent arrays in CSRMat constructor");
+for(SPMatindex i=0; i<n; ++i)
+	{
+	if(rp[i+1]<rp[i]) laerror("decreasing rowptr in CSRMat constructor");
+	for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) if(c[k]>=m || (k>rp[i] && c[k]<=c[k-1])) laerror("column index out of range or columns not ascending in CSRMat constructor");
+	}
+#endif
+}
+
+template <typename T>
+void CSRMat<T>::copyonwrite(bool detachonly, bool deep)
+{
+if(rowptr.size()) rowptr.copyonwrite(detachonly);
+if(col.size()) col.copyonwrite(detachonly);
+if(v.size()) v.copyonwrite(detachonly,deep);
+}
+
+template <typename T>
+void CSRMat<T>::resize(const SPMatindex n, const SPMatindex m)
+{
+nn=n; mm=m; symmetric=false;
+col.resize(0); v.resize(0);
+if(n==0 && m==0) {rowptr.resize(0); return;}
+rowptr.resize(n+1);
+rowptr.clear();
+}
+
+//counting sort, O(nnz+n+m); the columns of the result are ascending in each row also if they were not in the input
+template <typename T>
+const CSRMat<T> CSRMat<T>::transpose(bool conj) const
+{
+CSRMat<T> r(mm,nn);
+r.symmetric=symmetric;
+const int nz=v.size();
+if(nz==0) return r;
+const SPMatindex *rp=s_ptr(rowptr), *c=s_ptr(col);
+const T *x=s_ptr(v);
+r.col.resize(nz); r.v.resize(nz);
+SPMatindex *rrp=r.rowptr, *rc=r.col;
+T *rv=r.v;
+for(int k=0; k<nz; ++k) ++rrp[c[k]+1];
+for(SPMatindex j=0; j<mm; ++j) rrp[j+1]+=rrp[j];
+for(SPMatindex i=0; i<nn; ++i) for(SPMatindex k=rp[i]; k<rp[i+1]; ++k)
+	{
+	SPMatindex p=rrp[c[k]]++;
+	rc[p]=i;
+	rv[p]= conj? LA_traits<T>::conjugate(x[k]) : x[k];
+	}
+for(SPMatindex j=mm; j>0; --j) rrp[j]=rrp[j-1]; //rrp[j] was advanced to the start of row j+1
+rrp[0]=0;
+return r;
+}
+
+//removes small elements; also sums duplicates (adjacent in sorted rows, only during the construction from coordinates)
+template <typename T>
+void CSRMat<T>::simplify(const typename LA_traits<T>::normtype eps)
+{
+const int nz=v.size();
+if(nz==0) return;
+copyonwrite();
+SPMatindex *rp=rowptr, *c=col;
+T *x=v;
+int p=0, k=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	const int end=rp[i+1];
+	while(k<end)
+		{
+		const SPMatindex cc=c[k];
+		T s=x[k++];
+		while(k<end && c[k]==cc) s+=x[k++];
+		if(std::abs(s)>eps) {c[p]=cc; x[p]=s; ++p;}
+		}
+	rp[i+1]=p;
+	}
+if(p<nz) {col.resize(p,true); v.resize(p,true);}
+}
+
+//from coordinate format: grouped by columns (i.e. the transposed matrix in CSR, rows unsorted), transposed back (rows sorted), duplicates summed
+template <typename T>
+void CSRMat<T>::fromcoo(const SPMatindex n, const SPMatindex m, const int nz, const SPMatindex *ri, const SPMatindex *ci, const T *vi)
+{
+CSRMat<T> t(m,n);
+if(nz)
+	{
+	t.col.resize(nz); t.v.resize(nz);
+	SPMatindex *trp=t.rowptr, *tc=t.col;
+	T *tv=t.v;
+	for(int k=0; k<nz; ++k)
+		{
+#ifdef DEBUG
+		if(ri[k]>=n || ci[k]>=m) laerror("index out of range in CSRMat construction");
+#endif
+		++trp[ci[k]+1];
+		}
+	for(SPMatindex j=0; j<m; ++j) trp[j+1]+=trp[j];
+	for(int k=0; k<nz; ++k) {SPMatindex p=trp[ci[k]]++; tc[p]=ri[k]; tv[p]=vi[k];}
+	for(SPMatindex j=m; j>0; --j) trp[j]=trp[j-1];
+	trp[0]=0;
+	}
+*this=t.transpose();
+simplify();
+}
+
+//duplicates summed; a symmetric SparseMat (one triangle, complex symmetric) is unfolded without conjugation
+template <typename T>
+CSRMat<T>::CSRMat(const SparseMat<T> &rhs)
+: nn(0), mm(0), symmetric(false)
+{
+const bool sym=rhs.issymmetric();
+int nz=0;
+for(const matel<T> *l=rhs.getlist(); l; l=l->next) nz+= (sym && l->row!=l->col)? 2:1;
+NRVec<SPMatindex> ri(nz), ci(nz);
+NRVec<T> vi(nz);
+SPMatindex *r=s_wptr(ri), *c=s_wptr(ci);
+T *x=s_wptr(vi);
+int p=0;
+for(const matel<T> *l=rhs.getlist(); l; l=l->next)
+	{
+	r[p]=l->row; c[p]=l->col; x[p++]=l->elem;
+	if(sym && l->row!=l->col) {r[p]=l->col; c[p]=l->row; x[p++]=l->elem;}
+	}
+fromcoo(rhs.nrows(),rhs.ncols(),nz,r,c,x);
+symmetric= sym && !LA_traits<T>::is_complex();
+}
+
+template <typename T>
+CSRMat<T>::CSRMat(const SparseSMat<T> &rhs)
+: nn(rhs.nrows()), mm(rhs.ncols()), symmetric(false), rowptr(nn+1)
+{
+const typename LA_traits<T>::normtype eps=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon();
+typedef typename std::map<SPMatindex,T>::const_iterator IT;
+int nz=0;
+for(SPMatindex i=0; i<nn; ++i) {const std::map<SPMatindex,T> *l=rhs.line(i); if(l) for(IT q=l->begin(); q!=l->end(); ++q) if(std::abs(q->second)>eps) ++nz;}
+col.resize(nz); v.resize(nz);
+SPMatindex *rp=rowptr, *c=s_wptr(col);
+T *x=s_wptr(v);
+int p=0;
+rp[0]=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	const std::map<SPMatindex,T> *l=rhs.line(i);
+	if(l) for(IT q=l->begin(); q!=l->end(); ++q) if(std::abs(q->second)>eps) {c[p]=q->first; x[p++]=q->second;}
+	rp[i+1]=p;
+	}
+}
+
+template <typename T>
+CSRMat<T>::CSRMat(const NRMat<T> &rhs)
+: nn(rhs.nrows()), mm(rhs.ncols()), symmetric(false), rowptr(nn+1)
+{
+const typename LA_traits<T>::normtype eps=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon();
+int nz=0;
+for(SPMatindex i=0; i<nn; ++i) {const T *a=rhs[i]; for(SPMatindex j=0; j<mm; ++j) if(std::abs(a[j])>eps) ++nz;}
+col.resize(nz); v.resize(nz);
+SPMatindex *rp=rowptr, *c=s_wptr(col);
+T *x=s_wptr(v);
+int p=0;
+rp[0]=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	const T *a=rhs[i];
+	for(SPMatindex j=0; j<mm; ++j) if(std::abs(a[j])>eps) {c[p]=j; x[p++]=a[j];}
+	rp[i+1]=p;
+	}
+}
+
+//hermitian as NRSMat: the upper triangle is the conjugate of the stored lower one
+template <typename T>
+CSRMat<T>::CSRMat(const NRSMat<T> &rhs)
+: nn(rhs.nrows()), mm(rhs.ncols()), symmetric(true), rowptr(nn+1)
+{
+const typename LA_traits<T>::normtype eps=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon();
+int nz=0;
+for(SPMatindex i=0; i<nn; ++i) for(SPMatindex j=0; j<nn; ++j) if(std::abs(rhs(i,j))>eps) ++nz;
+col.resize(nz); v.resize(nz);
+SPMatindex *rp=rowptr, *c=s_wptr(col);
+T *x=s_wptr(v);
+int p=0;
+rp[0]=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	for(SPMatindex j=0; j<nn; ++j)
+		{
+		const T a=rhs(i,j);
+		if(std::abs(a)>eps) {c[p]=j; x[p++]= j>i? LA_traits<T>::conjugate(a) : a;}
+		}
+	rp[i+1]=p;
+	}
+}
+
+template <typename T>
+const T CSRMat<T>::operator()(const SPMatindex i, const SPMatindex j) const
+{
+#ifdef DEBUG
+if(i>=nn || j>=mm) laerror("index out of range in CSRMat::operator()");
+#endif
+if(v.size()==0) return (T)0;
+const SPMatindex *rp=rowptr, *c=col;
+const T *x=v;
+const SPMatindex *q=std::lower_bound(c+rp[i],c+rp[i+1],j);
+if(q!=c+rp[i+1] && *q==j) return x[q-c];
+return (T)0;
+}
+
+template <typename T>
+void CSRMat<T>::put(int fd, bool dimensions, bool transp, bool orcaformat) const
+{
+if(orcaformat) laerror("ORCA format not supported for CSRMat");
+if(transp) {transpose().put(fd,dimensions,false); return;}
+errno=0;
+if(dimensions)
+	{
+	if(sizeof(SPMatindex)!=write(fd,&nn,sizeof(SPMatindex))) laerror("cannot write() 1 in CSRMat::put");
+	if(sizeof(SPMatindex)!=write(fd,&mm,sizeof(SPMatindex))) laerror("cannot write() 2 in CSRMat::put");
+	}
+int sym=symmetric;
+if(sizeof(int)!=write(fd,&sym,sizeof(int))) laerror("cannot write() 3 in CSRMat::put");
+rowptr.put(fd,true);
+col.put(fd,true);
+v.put(fd,true);
+}
+
+template <typename T>
+void CSRMat<T>::get(int fd, bool dimensions, bool transp, bool orcaformat)
+{
+if(orcaformat) laerror("ORCA format not supported for CSRMat");
+errno=0;
+SPMatindex dim[2];
+if(dimensions) {if(2*sizeof(SPMatindex)!=read(fd,&dim,2*sizeof(SPMatindex))) laerror("read() error 1 in CSRMat::get");}
+else {dim[0]= transp?mm:nn; dim[1]= transp?nn:mm;}
+int sym;
+if(sizeof(int)!=read(fd,&sym,sizeof(int))) laerror("read() error 2 in CSRMat::get");
+rowptr.get(fd,true);
+col.get(fd,true);
+v.get(fd,true);
+nn=dim[0]; mm=dim[1]; symmetric=sym;
+if(rowptr.size()!=(nn||mm ? (int)nn+1 : 0) || col.size()!=v.size()) laerror("inconsistent data in CSRMat::get");
+if(transp) transposeme();
+}
+
+template <typename T>
+CSRMat<T> & CSRMat<T>::operator*=(const T &a)
+{
+if(a==(T)0) {clear(); return *this;}
+if(LA_traits<T>::imagpart(a)!=0) symmetric=false;
+if(v.size()) v*=a;
+return *this;
+}
+
+template <typename T>
+CSRMat<T> & CSRMat<T>::operator=(const T &a)
+{
+if(nn!=mm) laerror("assignment of scalar to non-square CSRMat");
+resize(nn,mm);
+symmetric= LA_traits<T>::imagpart(a)==0;
+if(nn==0 || std::abs(a)<=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()) return *this;
+col.resize(nn); v.resize(nn);
+SPMatindex *rp=rowptr, *c=col;
+T *x=v;
+for(SPMatindex i=0; i<nn; ++i) {rp[i]=i; c[i]=i; x[i]=a;}
+rp[nn]=nn;
+return *this;
+}
+
+//in place if all diagonal elements are stored (only the values are copied if shared), otherwise via axpy
+template <typename T>
+CSRMat<T> & CSRMat<T>::operator+=(const T &a)
+{
+if(nn!=mm) laerror("addition of scalar to non-square CSRMat");
+if(a==(T)0) return *this;
+const SPMatindex *rp=s_ptr(rowptr), *c=s_ptr(col);
+bool all= v.size()>0;
+if(all) for(SPMatindex i=0; i<nn; ++i)
+	{
+	const SPMatindex *q=std::lower_bound(c+rp[i],c+rp[i+1],i);
+	if(q==c+rp[i+1] || *q!=i) {all=false; break;}
+	}
+if(!all)
+	{
+	CSRMat<T> d(nn);
+	d=a;
+	const bool sym=symmetric;
+	axpy((T)1,d);
+	symmetric= sym && LA_traits<T>::imagpart(a)==0;
+	return *this;
+	}
+if(LA_traits<T>::imagpart(a)!=0) symmetric=false;
+v.copyonwrite();
+T *x=v;
+for(SPMatindex i=0; i<nn; ++i) x[std::lower_bound(c+rp[i],c+rp[i+1],i)-c]+=a;
+return *this;
+}
+
+//merge of the sorted rows; for the same structure (shared or equal arrays) only the values by BLAS axpy
+//the result is declared hermitian if both are (a matrix without elements counts as hermitian) and alpha is real
+template <typename T>
+void CSRMat<T>::axpy(const T alpha, const CSRMat<T> &x0, const bool transp)
+{
+if((!transp && (nn!=x0.nn || mm!=x0.mm)) || (transp && (nn!=x0.mm || mm!=x0.nn))) laerror("incompatible dimensions for axpy");
+if(alpha==(T)0 || x0.v.size()==0) return;
+CSRMat<T> xt;
+if(transp) xt=x0.transpose();
+const CSRMat<T> &x= transp? xt : x0;
+const bool sym= (v.size()==0? nn==mm : symmetric) && x.symmetric && LA_traits<T>::imagpart(alpha)==0;
+if(v.size()==0)
+	{
+	rowptr=x.rowptr; col=x.col; v=x.v;
+	if(alpha!=(T)1) v*=alpha;
+	symmetric=sym;
+	return;
+	}
+const SPMatindex *rp=rowptr, *c=col, *xrp=x.rowptr, *xc=x.col;
+const int nz=v.size(), xnz=x.v.size();
+if(nz==xnz && ((rp==xrp && c==xc) || (!memcmp(rp,xrp,(nn+1)*sizeof(SPMatindex)) && !memcmp(c,xc,nz*sizeof(SPMatindex)))))
+	{
+	v.axpy(alpha,x.v);
+	symmetric=sym;
+	return;
+	}
+const T *a=v, *b=x.v;
+NRVec<SPMatindex> nrp(nn+1);
+SPMatindex *np=nrp;
+np[0]=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	SPMatindex p=rp[i], q=xrp[i], cnt=0;
+	while(p<rp[i+1] && q<xrp[i+1]) {if(c[p]<xc[q]) ++p; else if(c[p]>xc[q]) ++q; else {++p; ++q;} ++cnt;}
+	np[i+1]=np[i]+cnt+(rp[i+1]-p)+(xrp[i+1]-q);
+	}
+NRVec<SPMatindex> ncol(np[nn]);
+NRVec<T> nv(np[nn]);
+SPMatindex *nc=ncol;
+T *y=nv;
+int k=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	SPMatindex p=rp[i], q=xrp[i];
+	while(p<rp[i+1] || q<xrp[i+1])
+		{
+		if(q==xrp[i+1] || (p<rp[i+1] && c[p]<xc[q])) {nc[k]=c[p]; y[k++]=a[p++];}
+		else if(p==rp[i+1] || c[p]>xc[q]) {nc[k]=xc[q]; y[k++]=alpha*b[q++];}
+		else {nc[k]=c[p]; y[k++]=a[p++]+alpha*b[q++];}
+		}
+	}
+rowptr=nrp; col=ncol; v=nv;
+symmetric=sym;
+}
+
+template <typename T>
+const typename LA_traits<T>::normtype CSRMat<T>::norm(const T scalar) const
+{
+if(scalar==(T)0) return v.size()? v.norm() : (typename LA_traits<T>::normtype)0;
+if(nn!=mm) laerror("subtraction of scalar from non-square CSRMat in norm()");
+const SPMatindex *rp=s_ptr(rowptr), *c=s_ptr(col);
+const T *x=s_ptr(v);
+typename LA_traits<T>::normtype sum=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	bool d=false;
+	if(x) for(SPMatindex k=rp[i]; k<rp[i+1]; ++k)
+		{
+		T e=x[k];
+		if(c[k]==i) {e-=scalar; d=true;}
+		sum+=LA_traits<T>::sqrabs(e);
+		}
+	if(!d) sum+=LA_traits<T>::sqrabs(scalar);
+	}
+return std::sqrt(sum);
+}
+
+template <typename T>
+const T CSRMat<T>::dot(const CSRMat<T> &rhs) const
+{
+if(nn!=rhs.nn || mm!=rhs.mm) laerror("dot of incompatible CSRMat matrices");
+if(v.size()==0 || rhs.v.size()==0) return (T)0;
+const SPMatindex *rp=rowptr, *c=col, *xrp=rhs.rowptr, *xc=rhs.col;
+if(rp==xrp && c==xc) return v.dot(rhs.v);
+const T *a=v, *b=rhs.v;
+T s=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	SPMatindex p=rp[i], q=xrp[i];
+	while(p<rp[i+1] && q<xrp[i+1])
+		{
+		if(c[p]<xc[q]) ++p;
+		else if(c[p]>xc[q]) ++q;
+		else s+=LA_traits<T>::conjugate(a[p++])*b[q++];
+		}
+	}
+return s;
+}
+
+template <typename T>
+const T CSRMat<T>::dot(const NRMat<T> &rhs) const
+{
+if((int)nn!=rhs.nrows() || (int)mm!=rhs.ncols()) laerror("dot of incompatible matrices");
+if(v.size()==0) return (T)0;
+const SPMatindex *rp=rowptr, *c=col;
+const T *a=v;
+T s=0;
+for(SPMatindex i=0; i<nn; ++i)
+	{
+	const T *b=rhs[i];
+	for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) s+=LA_traits<T>::conjugate(a[k])*b[c[k]];
+	}
+return s;
+}
+
+template <typename T>
+const T CSRMat<T>::trace() const
+{
+if(v.size()==0) return (T)0;
+const SPMatindex *rp=rowptr, *c=col;
+const T *x=v;
+T s=0;
+for(SPMatindex i=0; i<nn && i<mm; ++i)
+	{
+	const SPMatindex *q=std::lower_bound(c+rp[i],c+rp[i+1],i);
+	if(q!=c+rp[i+1] && *q==i) s+=x[q-c];
+	}
+return s;
+}
+
+//sparse times sparse, row-wise with a dense accumulator: F. G. Gustavson, ACM Trans. Math. Softw. 4, 250 (1978)
+//two passes (count, then compute); the columns of each result row are sorted; the result is not declared hermitian
+template <typename T>
+void CSRMat<T>::gemm(const T beta, const CSRMat<T> &a0, const char transa, const CSRMat<T> &b0, const char transb, const T alpha)
+{
+const char ta=tolower(transa), tb=tolower(transb);
+CSRMat<T> at, bt;
+if(ta!='n') at=a0.transpose(ta=='c');
+if(tb!='n') bt=b0.transpose(tb=='c');
+const CSRMat<T> &a= ta!='n'? at : a0;
+const CSRMat<T> &b= tb!='n'? bt : b0;
+if(a.mm!=b.nn) laerror("incompatible matrices in CSRMat::gemm");
+if(beta!=(T)0 && (nn!=a.nn || mm!=b.mm)) laerror("incompatible result matrix in CSRMat::gemm");
+const SPMatindex n=a.nn, m=b.mm;
+CSRMat<T> r(n,m);
+if(alpha!=(T)0 && a.v.size() && b.v.size())
+	{
+	const SPMatindex *arp=a.rowptr, *ac=a.col, *brp=b.rowptr, *bc=b.col;
+	const T *av=a.v, *bv=b.v;
+	NRVec<SPMatindex> marker(m);
+	SPMatindex *mk=marker;
+	for(SPMatindex j=0; j<m; ++j) mk[j]=n;
+	SPMatindex *rp=r.rowptr;
+	for(SPMatindex i=0; i<n; ++i)
+		{
+		SPMatindex cnt=0;
+		for(SPMatindex k=arp[i]; k<arp[i+1]; ++k) for(SPMatindex q=brp[ac[k]]; q<brp[ac[k]+1]; ++q) if(mk[bc[q]]!=i) {mk[bc[q]]=i; ++cnt;}
+		rp[i+1]=rp[i]+cnt;
+		}
+	r.col.resize(rp[n]); r.v.resize(rp[n]);
+	if(rp[n])
+		{
+		SPMatindex *c=r.col;
+		T *y=r.v;
+		NRVec<T> acc(m);
+		T *ac2=acc;
+		for(SPMatindex j=0; j<m; ++j) mk[j]=n;
+		for(SPMatindex i=0; i<n; ++i)
+			{
+			SPMatindex len=0;
+			for(SPMatindex k=arp[i]; k<arp[i+1]; ++k)
+				{
+				const T x=av[k];
+				for(SPMatindex q=brp[ac[k]]; q<brp[ac[k]+1]; ++q)
+					{
+					const SPMatindex j=bc[q];
+					if(mk[j]!=i) {mk[j]=i; c[rp[i]+len++]=j; ac2[j]=x*bv[q];}
+					else ac2[j]+=x*bv[q];
+					}
+				}
+			std::sort(c+rp[i],c+rp[i]+len);
+			for(SPMatindex t=rp[i]; t<rp[i]+len; ++t) y[t]=alpha*ac2[c[t]];
+			}
+		}
+	}
+if(beta==(T)0) *this=r;
+else
+	{
+	*this*=beta;
+	axpy((T)1,r);
+	}
+symmetric=false;
+}
+
+template <typename T>
+CSRMat<T> & CSRMat<T>::oplusequal(const CSRMat<T> &rhs)
+{
+const bool sym= (v.size()==0? nn==mm : symmetric) && rhs.symmetric;
+const SPMatindex n=nn+rhs.nn, m=mm+rhs.mm;
+const int nz=v.size(), rnz=rhs.v.size();
+NRVec<SPMatindex> nrp(n+1), ncol(nz+rnz);
+NRVec<T> nv(nz+rnz);
+SPMatindex *np=nrp, *nc=s_wptr(ncol);
+T *y=s_wptr(nv);
+const SPMatindex *rp=s_ptr(rowptr), *c=s_ptr(col), *rrp=s_ptr(rhs.rowptr), *rc=s_ptr(rhs.col);
+const T *x=s_ptr(v), *rx=s_ptr(rhs.v);
+for(SPMatindex i=0; i<=nn; ++i) np[i]= rp? rp[i] : 0;
+for(SPMatindex i=1; i<=rhs.nn; ++i) np[nn+i]=nz+rrp[i];
+for(int k=0; k<nz; ++k) {nc[k]=c[k]; y[k]=x[k];}
+for(int k=0; k<rnz; ++k) {nc[nz+k]=mm+rc[k]; y[nz+k]=rx[k];}
+nn=n; mm=m;
+rowptr=nrp; col=ncol; v=nv;
+symmetric=sym;
+return *this;
+}
+
+template <typename T>
+const CSRMat<T> CSRMat<T>::otimes(const CSRMat<T> &rhs) const
+{
+CSRMat<T> r(nn*rhs.nn,mm*rhs.mm);
+r.symmetric= symmetric && rhs.symmetric;
+const int nz=v.size(), rnz=rhs.v.size();
+if(nz==0 || rnz==0) return r;
+const SPMatindex *rp=rowptr, *c=col, *brp=rhs.rowptr, *bc=rhs.col;
+const T *x=v, *bx=rhs.v;
+r.col.resize(nz*rnz); r.v.resize(nz*rnz);
+SPMatindex *np=r.rowptr, *nc=r.col;
+T *y=r.v;
+int p=0;
+for(SPMatindex i=0; i<nn; ++i) for(SPMatindex k=0; k<rhs.nn; ++k)
+	{
+	for(SPMatindex a=rp[i]; a<rp[i+1]; ++a) for(SPMatindex b=brp[k]; b<brp[k+1]; ++b) {nc[p]=c[a]*rhs.mm+bc[b]; y[p++]=x[a]*bx[b];}
+	np[i*rhs.nn+k+1]=p;
+	}
+return r;
+}
+
+//row r becomes row to[r]: the rows are moved as blocks, their columns stay sorted
+template <typename T>
+void CSRMat<T>::permrows(const SPMatindex *to)
+{
+const int nz=v.size();
+NRVec<SPMatindex> nrp(nn+1);
+SPMatindex *np=nrp;
+const SPMatindex *rp=rowptr;
+np[0]=0;
+for(SPMatindex r=0; r<nn; ++r) np[to[r]+1]=rp[r+1]-rp[r];
+for(SPMatindex r=0; r<nn; ++r) np[r+1]+=np[r];
+if(nz)
+	{
+	NRVec<SPMatindex> ncol(nz);
+	NRVec<T> nv(nz);
+	SPMatindex *nc=ncol;
+	T *y=nv;
+	const SPMatindex *c=col;
+	const T *x=v;
+	for(SPMatindex r=0; r<nn; ++r) for(SPMatindex k=rp[r], t=np[to[r]]; k<rp[r+1]; ++k, ++t) {nc[t]=c[k]; y[t]=x[k];}
+	col=ncol; v=nv;
+	}
+rowptr=nrp;
+}
+
+template <typename T>
+void CSRMat<T>::permuteme_rows(const NRPerm<int> &p, const bool inverse)
+{
+if((SPMatindex)p.size()!=nn) laerror("inconsistent dimension in permuteme_rows");
+NRPerm<int> pinv;
+if(inverse) pinv=p.inverse();
+const NRPerm<int> &pp= inverse?pinv:p;
+NRVec<SPMatindex> to(nn);
+for(SPMatindex r=0; r<nn; ++r) to[r]=pp[r+1]-1;
+if(nn) permrows(to);
+symmetric=false;
+}
+
+//columns via the transposed matrix (counting sorts keep the rows sorted)
+template <typename T>
+void CSRMat<T>::permuteme_cols(const NRPerm<int> &p, const bool inverse)
+{
+if((SPMatindex)p.size()!=mm) laerror("inconsistent dimension in permuteme_cols");
+*this=transpose();
+permuteme_rows(p,inverse);
+*this=transpose();
+}
+
+template <typename T>
+void CSRMat<T>::permuteme_both(const NRPerm<int> &p, const NRPerm<int> &q, const bool inverse)
+{
+const bool sym= symmetric && p==q;
+permuteme_rows(p,inverse);
+permuteme_cols(q,inverse);
+symmetric=sym;
+}
+
+template <typename T>
+void CSRMat<T>::permuterows(const NRVec<SPMatindex> &p)
+{
+if((SPMatindex)p.size()!=nn) laerror("inconsistent dimension in permuterows");
+if(nn) permrows(p);
+symmetric=false;
+}
+
+template <typename T>
+void CSRMat<T>::permutecolumns(const NRVec<SPMatindex> &p)
+{
+if((SPMatindex)p.size()!=mm) laerror("inconsistent dimension in permutecolumns");
+*this=transpose();
+permuterows(p);
+*this=transpose();
+}
+
+template <typename T>
+void CSRMat<T>::permuteindices(const NRVec<SPMatindex> &p)
+{
+if((SPMatindex)p.size()!=nn || (SPMatindex)p.size()!=mm) laerror("inconsistent dimension in permuteindices");
+const bool sym=symmetric;
+permuterows(p);
+permutecolumns(p);
+symmetric=sym;
+}
+
+//square: diagonal (binary search in each row); non-square: diagonal of A^T A (for preconditioning as SparseMat::diagonalof); cache is ignored
+template <typename T>
+const T* CSRMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const
+{
+#ifdef DEBUG
+if((int)mm!=r.size()) laerror("incompatible vector size in diagonalof()");
+#endif
+NOT_GPU(r);
+r.copyonwrite();
+T *rr=r;
+const SPMatindex *rp=s_ptr(rowptr), *c=s_ptr(col);
+const T *x=s_ptr(v);
+if(nn==mm)
+	{
+	for(SPMatindex i=0; i<nn; ++i)
+		{
+		T d=0;
+		if(x)
+			{
+			const SPMatindex *q=std::lower_bound(c+rp[i],c+rp[i+1],i);
+			if(q!=c+rp[i+1] && *q==i) d=x[q-c];
+			}
+		if(divide) {if(d!=(T)0) rr[i]/=d;}
+		else rr[i]=d;
+		}
+	}
+else
+	{
+	NRVec<T> d(mm);
+	d.clear();
+	T *dd=d;
+	for(int k=0; k<v.size(); ++k) dd[c[k]]+=x[k]*x[k];
+	if(divide) {for(SPMatindex j=0; j<mm; ++j) if(dd[j]!=(T)0) rr[j]/=dd[j];}
+	else for(SPMatindex j=0; j<mm; ++j) rr[j]=dd[j];
+	}
+return divide?NULL:rr;
+}
+
+//r = beta r + alpha op(A) x, op = 'n', 't' or 'c' (conjugate transpose, the same as 't' for real types)
+template <typename T>
+void NRVec<T>::gemv(const T beta, const CSRMat<T> &a, const char trans, const T alpha, const NRVec<T> &x)
+{
+const char tr=tolower(trans);
+const bool tp= tr!='n';
+const bool cj= tr=='c' && LA_traits<T>::is_complex();
+if((tp?a.nrows():a.ncols())!=(SPMatindex)x.size()) laerror("incompatible sizes in gemv");
+#ifdef DEBUG
+if((tp?a.ncols():a.nrows())!=(SPMatindex)nn) laerror("incompatible result vector in gemv");
+#endif
+NOT_GPU(*this);
+NOT_GPU(x);
+copyonwrite();
+if(beta==(T)0) clear();
+else if(beta!=(T)1) (*this)*=beta;
+if(alpha==(T)0 || a.length()==0) return;
+const SPMatindex *rp=a.getrowptr(), *c=a.getcol();
+const T *av=a.getval();
+const T *xv=x.v;
+const SPMatindex n=a.nrows();
+if(!tp)
+	for(SPMatindex i=0; i<n; ++i)
+		{
+		T s=0;
+		for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) s+=av[k]*xv[c[k]];
+		v[i]+=alpha*s;
+		}
+else
+	for(SPMatindex i=0; i<n; ++i)
+		{
+		const T ax=alpha*xv[i];
+		if(cj) for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) v[c[k]]+=LA_traits<T>::conjugate(av[k])*ax;
+		else for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) v[c[k]]+=av[k]*ax;
+		}
+}
+
+//vector from a matrix with one row or one column
+template <typename T>
+NRVec<T>::NRVec(const CSRMat<T> &rhs)
+: NRVec(rhs.nrows()>1?rhs.nrows():rhs.ncols())
+{
+if(rhs.nrows()>1 && rhs.ncols()>1) laerror("cannot construct a vector from a sparse matrix with more than one row/column");
+clear();
+const SPMatindex *rp=s_ptr(rhs.getrowptr()), *c=s_ptr(rhs.getcol());
+const T *x=s_ptr(rhs.getval());
+if(!x) return;
+if(rhs.nrows()>1) {for(SPMatindex i=0; i<rhs.nrows(); ++i) if(rp[i+1]>rp[i]) v[i]=x[rp[i]];}
+else for(SPMatindex k=0; k<rp[1]; ++k) v[c[k]]=x[k];
+}
+
+//conversions to the other classes
+
+template <typename T>
+NRMat<T>::NRMat(const CSRMat<T> &rhs)
+{
+nn=rhs.nrows();
+mm=rhs.ncols();
+count=new int(1);
+#ifdef CUDALA
+location=cpu;
+#endif
+T *p;
+#ifdef MATPTR
+	v= new T*[nn];
+	p=v[0] = new T[mm*nn];
+	for (int i=1; i< nn; i++) v[i] = v[i-1] + mm;
+#else
+	p= v = new T[mm*nn];
+#endif
+memset(p,0,nn*mm*sizeof(T));
+const SPMatindex *rp=s_ptr(rhs.getrowptr()), *c=s_ptr(rhs.getcol());
+const T *x=s_ptr(rhs.getval());
+if(x) for(int i=0; i<nn; ++i) for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) p[(size_t)i*mm+c[k]]=x[k];
+}
+
+template <typename T>
+NRMat<T> & NRMat<T>::operator+=(const CSRMat<T> &rhs)
+{
+if(nn!=(int)rhs.nrows() || mm!=(int)rhs.ncols()) laerror("incompatible matrices in +=");
+NOT_GPU(*this);
+copyonwrite();
+const SPMatindex *rp=s_ptr(rhs.getrowptr()), *c=s_ptr(rhs.getcol());
+const T *x=s_ptr(rhs.getval());
+if(x) for(int i=0; i<nn; ++i) {T *a=(*this)[i]; for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) a[c[k]]+=x[k];}
+return *this;
+}
+
+template <typename T>
+NRMat<T> & NRMat<T>::operator-=(const CSRMat<T> &rhs)
+{
+if(nn!=(int)rhs.nrows() || mm!=(int)rhs.ncols()) laerror("incompatible matrices in -=");
+NOT_GPU(*this);
+copyonwrite();
+const SPMatindex *rp=s_ptr(rhs.getrowptr()), *c=s_ptr(rhs.getcol());
+const T *x=s_ptr(rhs.getval());
+if(x) for(int i=0; i<nn; ++i) {T *a=(*this)[i]; for(SPMatindex k=rp[i]; k<rp[i+1]; ++k) a[c[k]]-=x[k];}
+return *this;
+}
+
+//C = beta C + alpha op(A) op(B), A sparse: rows of op(B) accumulated (op(B) transposed first if needed)
+template <typename T>
+void NRMat<T>::gemm(const T &beta, const CSRMat<T> &a, const char transa, const NRMat<T> &b, const char transb, const T &alpha)
+{
+const char ta=tolower(transa), tb=tolower(transb);
+const bool tpa= ta!='n', cja= ta=='c' && LA_traits<T>::is_complex();
+const int an= tpa? a.ncols() : a.nrows(), ak= tpa? a.nrows() : a.ncols();
+const int bk= tb!='n'? b.ncols() : b.nrows(), bm= tb!='n'? b.nrows() : b.ncols();
+if(ak!=bk) laerror("incompatible matrices in NRMat::gemm with CSRMat");
+#ifdef DEBUG
+if(nn!=an || mm!=bm) laerror("incompatible result matrix in NRMat::gemm with CSRMat");
+#endif
+NOT_GPU(*this);
+copyonwrite();
+if(beta==(T)0) clear();
+else if(beta!=(T)1) (*this)*=beta;
+if(alpha==(T)0 || a.length()==0) return;
+NRMat<T> bt;
+if(tb!='n') bt=b.transpose(tb=='c');
+const NRMat<T> &bb= tb!='n'? bt : b;
+const SPMatindex *rp=a.getrowptr(), *c=a.getcol();
+const T *x=a.getval();
+for(SPMatindex i=0; i<a.nrows(); ++i) for(SPMatindex k=rp[i]; k<rp[i+1]; ++k)
+	{
+	const T s= alpha*(cja? LA_traits<T>::conjugate(x[k]) : x[k]);
+	T *d= (*this)[tpa? c[k] : i];
+	const T *e= bb[tpa? i : c[k]];
+	for(int j=0; j<mm; ++j) d[j]+=s*e[j];
+	}
+}
+
+template <typename T>
+const NRMat<T> NRMat<T>::operator*(const CSRMat<T> &rhs) const
+{
+if((SPMatindex)mm!=rhs.nrows()) laerror("incompatible matrices in NRMat*CSRMat");
+NOT_GPU(*this);
+NRMat<T> r(nn,rhs.ncols());
+r.clear();
+if(rhs.length()==0) return r;
+const SPMatindex *rp=rhs.getrowptr(), *c=rhs.getcol();
+const T *x=rhs.getval();
+for(int i=0; i<nn; ++i)
+	{
+	const T *a=(*this)[i];
+	T *d=r[i];
+	for(int k=0; k<mm; ++k) if(a[k]!=(T)0) for(SPMatindex q=rp[k]; q<rp[k+1]; ++q) d[c[q]]+=a[k]*x[q];
+	}
+return r;
+}
+
+//lower triangle, as from SparseSMat
+template <typename T>
+NRSMat<T>::NRSMat(const CSRMat<T> &rhs)
+: nn(rhs.nrows())
+{
+if(rhs.nrows()!=rhs.ncols()) laerror("cannot transform rectangular matrix to NRSMat");
+#ifdef CUDALA
+location=cpu;
+#endif
+count=new int(1);
+v=new T[NN2];
+memset(v,0,NN2*sizeof(T));
+const SPMatindex *rp=s_ptr(rhs.getrowptr()), *c=s_ptr(rhs.getcol());
+const T *x=s_ptr(rhs.getval());
+if(x) for(int i=0; i<nn; ++i) for(SPMatindex k=rp[i]; k<rp[i+1] && (int)c[k]<=i; ++k) v[(size_t)i*(i+1)/2+c[k]]=x[k];
+}
+
+template <typename T>
+SparseMat<T>::SparseMat(const CSRMat<T> &rhs)
+: nn(rhs.nrows()), mm(rhs.ncols()), symmetric(0), nonzero(0), count(new int(1)), list(NULL), rowsorted(NULL), colsorted(NULL)
+{
+for(typename CSRMat<T>::iterator p=rhs.begin(); p.notend(); ++p) add(p->row,p->col,p->elem);
+}
+
+template <typename T>
+SparseSMat<T>::SparseSMat(const CSRMat<T> &rhs)
+: nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1))
+{
+v= new std::map<SPMatindex,T> * [nn];
+memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
+for(typename CSRMat<T>::iterator p=rhs.begin(); p.notend(); ++p) add(p->row,p->col,p->elem,false);
+}


+//real and imaginary parts, complex matrices from real ones (nonclass.h): the index arrays stay shared,
+//unless elements which became zero have to be removed
+template<typename R>
+static const CSRMat<R> s_part(const CSRMat<std::complex<R> > &a, const NRVec<R> &x, const bool sym)
+{
+CSRMat<R> r(a.nrows(),a.ncols());
+if(a.length())
+	{
+	r=CSRMat<R>(a.nrows(),a.ncols(),a.getrowptr(),a.getcol(),x);
+	const R eps=LA_precision<R>::sparse_epsilon();
+	const R *y=x;
+	for(int k=0; k<x.size(); ++k) if(std::abs(y[k])<=eps) {r.simplify(); break;}
+	}
+if(sym && a.nrows()==a.ncols()) r.setsymmetric();
+return r;
+}
+
+template<typename R>
+static const CSRMat<std::complex<R> > s_cmplx(const CSRMat<R> &a, const bool imag)
+{
+CSRMat<std::complex<R> > r(a.nrows(),a.ncols());
+if(a.length()) r=CSRMat<std::complex<R> >(a.nrows(),a.ncols(),a.getrowptr(),a.getcol(),NRVec<std::complex<R> >(a.getval(),imag));
+if(!imag && a.issymmetric()) r.setsymmetric();
+return r;
+}
+
+template<>
+const CSRMat<double> realpart<CSRMat<std::complex<double> > >(const CSRMat<std::complex<double> > &a)
+{return s_part(a,a.length()? realpart(a.getval()) : NRVec<double>(),a.issymmetric());}
+
+template<>
+const CSRMat<float> realpart<CSRMat<std::complex<float> > >(const CSRMat<std::complex<float> > &a)
+{return s_part(a,a.length()? realpart(a.getval()) : NRVec<float>(),a.issymmetric());}
+
+template<>
+const CSRMat<double> imagpart<CSRMat<std::complex<double> > >(const CSRMat<std::complex<double> > &a)
+{return s_part(a,a.length()? imagpart(a.getval()) : NRVec<double>(),false);}
+
+template<>
+const CSRMat<float> imagpart<CSRMat<std::complex<float> > >(const CSRMat<std::complex<float> > &a)
+{return s_part(a,a.length()? imagpart(a.getval()) : NRVec<float>(),false);}
+
+template<>
+const CSRMat<std::complex<double> > realmatrix<CSRMat<double> >(const CSRMat<double> &a)
+{return s_cmplx(a,false);}
+
+template<>
+const CSRMat<std::complex<float> > realmatrix<CSRMat<float> >(const CSRMat<float> &a)
+{return s_cmplx(a,false);}
+
+template<>
+const CSRMat<std::complex<double> > imagmatrix<CSRMat<double> >(const CSRMat<double> &a)
+{return s_cmplx(a,true);}
+
+template<>
+const CSRMat<std::complex<float> > imagmatrix<CSRMat<float> >(const CSRMat<float> &a)
+{return s_cmplx(a,true);}
+
+//re + i im: a BLAS axpy of the values if both have the same structure, otherwise a merge
+template<>
+const CSRMat<std::complex<double> > complexmatrix<CSRMat<double> >(const CSRMat<double> &re, const CSRMat<double> &im)
+{CSRMat<std::complex<double> > r=s_cmplx(re,false); r.axpy(std::complex<double>(0,1),s_cmplx(im,false)); return r;}
+
+template<>
+const CSRMat<std::complex<float> > complexmatrix<CSRMat<float> >(const CSRMat<float> &re, const CSRMat<float> &im)
+{CSRMat<std::complex<float> > r=s_cmplx(re,false); r.axpy(std::complex<float>(0,1),s_cmplx(im,false)); return r;}

-/*
- Commented out by Roman for ICC

 #define INSTANTIZE(T) \
-template void CSRMat<T>::gemm(const T beta, const CSRMat &a, const char transa, const CSRMat &b, const char transb, const T alpha); \
-template CSRMat<T> & CSRMat<T>::operator*=(const T &a); \
-template void CSRMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const; \
-template void CSRMat<T>::axpy(const T alpha, const CSRMat &x, const bool transp); \
-template CSRMat<T> & CSRMat<T>::operator=(const T &a); \
-template CSRMat<T> & CSRMat<T>::operator+=(const T &a); \
-template CSRMat<T> & CSRMat<T>::operator-=(const T &a); \
-template LA_traits<T>::normtype CSRMat<T>::norm(const T scalar) const; \
-template const T* CSRMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \
-template void CSRMat<T>::get(int fd, bool dimen, bool transp); \
-template void CSRMat<T>::put(int fd, bool dimen, bool transp) const; \
-
+template class CSRMat<T>; \
+template NRMat<T>::NRMat(const CSRMat<T> &rhs); \
+template NRSMat<T>::NRSMat(const CSRMat<T> &rhs); \
+template SparseMat<T>::SparseMat(const CSRMat<T> &rhs); \
+template SparseSMat<T>::SparseSMat(const CSRMat<T> &rhs); \
+template void NRVec<T>::gemv(const T beta, const CSRMat<T> &a, const char trans, const T alpha, const NRVec<T> &x); \
+template NRVec<T>::NRVec(const CSRMat<T> &rhs); \
+template NRMat<T> & NRMat<T>::operator+=(const CSRMat<T> &rhs); \
+template NRMat<T> & NRMat<T>::operator-=(const CSRMat<T> &rhs); \
+template void NRMat<T>::gemm(const T &beta, const CSRMat<T> &a, const char transa, const NRMat<T> &b, const char transb, const T &alpha); \
+template const NRMat<T> NRMat<T>::operator*(const CSRMat<T> &rhs) const; \

 INSTANTIZE(double)
+INSTANTIZE(float)
 INSTANTIZE(std::complex<double>)
-*/
-
-//// forced instantization of functions in the header in the corresponding object file
-//@@@@@template class CSRMat<double>;
-//@@@@template class CSRMat<complex<double> >;
-//@@@@ unfinished class commented out
-
-
+INSTANTIZE(std::complex<float>)

 }//namespace
@@ -1,6 +1,6 @@
 /*
    LA: linear algebra C++ interface library
-    Copyright (C) 2010 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
+    Copyright (C) 2010-2026 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -32,16 +32,13 @@
 #include "vec.h"
 #include "mat.h"
 #include "smat.h"
-#include "qsort.h"
-
-#ifdef CUDALA
-#include "cusparse.h"
-#include "cuda_runtime.h"
-#endif

 namespace LA {

 //compressed row format read-efficient matrix storage
+//row i has the elements v[k], column col[k], k=rowptr[i]..rowptr[i+1]-1, columns strictly ascending within a row
+//the full matrix is always stored; symmetric only declares A = A^dagger (hermitian, real symmetric) e.g. for davidson, lanczos
+//rowptr, col and v share data independently (copy-on-write), e.g. a scaled copy shares the structure with the original

 template <typename T>
 class CSRMat
@@ -49,92 +46,157 @@ class CSRMat
 protected:
 	SPMatindex nn;
 	SPMatindex mm;
-	SPMatindex nnz;
-	int *count;
-	SPMatindex *rowptr;
-	SPMatindex *col;
-	T *v;
-#ifdef CUDALA
-        GPUID location;
-	cusparseMatDescr_t descr;
-#endif
-
+	bool symmetric;
+	NRVec<SPMatindex> rowptr;
+	NRVec<SPMatindex> col;
+	NRVec<T> v;
+	void fromcoo(const SPMatindex n, const SPMatindex m, const int nz, const SPMatindex *ri, const SPMatindex *ci, const T *vi); //sorts, sums duplicates, drops small elements
+	void permrows(const SPMatindex *to); //row r becomes row to[r]
 public:
-	CSRMat() : nn(0), mm(0), nnz(0), count(NULL), rowptr(NULL), col(NULL), v(NULL) {};
-	CSRMat(const CSRMat &rhs);
-	explicit CSRMat(const SPMatindex n, const SPMatindex m); 
-        explicit CSRMat(const SPMatindex n);
+	//iterator over the stored elements in row-major order
+	class iterator {
+	private:
+		const CSRMat *a;
+		SPMatindex row;
+		int k;
+		matel<T> my;
+		void set() {while(row<a->nn && k>=(int)a->rowptr[row+1]) ++row; if(row<a->nn) {my.row=row; my.col=a->col[k]; my.elem=a->v[k]; my.next=NULL;}}
+	public:
+		iterator(): a(NULL), row(0), k(0) {};
+		iterator(const CSRMat &rhs): a(&rhs), row(0), k(0) {if(a->nn) set();}
+		iterator & operator++() {++k; set(); return *this;}
+		bool notend() const {return a && row<a->nn;}
+		bool operator==(const iterator &rhs) const {return notend()==rhs.notend() && (!notend() || k==rhs.k);}
+		bool operator!=(const iterator &rhs) const {return !(*this==rhs);}
+		const matel<T> & operator*() const {return my;}
+		const matel<T> * operator->() const {return &my;}
+	};
+	iterator begin() const {return iterator(*this);}
+	iterator end() const {return iterator();}
+
+	CSRMat() : nn(0), mm(0), symmetric(false) {};
+	explicit CSRMat(const SPMatindex n, const SPMatindex m); //zero matrix
+	explicit CSRMat(const SPMatindex n); //zero square matrix
+	CSRMat(const SPMatindex n, const SPMatindex m, const NRVec<SPMatindex> &rowptr, const NRVec<SPMatindex> &col, const NRVec<T> &v); //shares the arrays
 	explicit CSRMat(const SparseMat<T> &rhs);
 	explicit CSRMat(const SparseSMat<T> &rhs);
 	explicit CSRMat(const NRSMat<T> &rhs);
 	explicit CSRMat(const NRMat<T> &rhs);
-	CSRMat & operator=(const CSRMat &rhs);
-	void copyonwrite();
-        void resize(const SPMatindex nn, const SPMatindex mm);
+	void copyonwrite(bool detachonly=false, bool deep=true);
+	void resize(const SPMatindex n, const SPMatindex m); //zero matrix
 	void dealloc(void) {resize(0,0);}
-	void moveto(GPUID destination);
-        void clear();
-	~CSRMat();
-	inline int getcount() const {return count?*count:0;}
+	void clear() {bool sym=symmetric; resize(nn,mm); symmetric=sym;} //zero matrix of the same dimensions and symmetry
+	void swap(CSRMat &rhs)
+		{
+		SPMatindex tmpnn=nn; nn=rhs.nn; rhs.nn=tmpnn;
+		SPMatindex tmpmm=mm; mm=rhs.mm; rhs.mm=tmpmm;
+		bool tmpsymmetric=symmetric; symmetric=rhs.symmetric; rhs.symmetric=tmpsymmetric;
+		rowptr.swap(rhs.rowptr); col.swap(rhs.col); v.swap(rhs.v);
+		}
+	inline SPMatindex nrows() const {return nn;}
+	inline SPMatindex ncols() const {return mm;}
+	inline int length() const {return v.size();} //number of stored elements
+	inline bool issymmetric() const {return symmetric;}
+	inline void setsymmetric() {if(nn!=mm) laerror("non-square cannot be symmetric"); symmetric=true;} //only declares the matrix hermitian
+	inline void setunsymmetric() {symmetric=false;}
+	inline const NRVec<SPMatindex> & getrowptr() const {return rowptr;}
+	inline const NRVec<SPMatindex> & getcol() const {return col;}
+	inline const NRVec<T> & getval() const {return v;}
+	const T operator()(const SPMatindex i, const SPMatindex j) const; //element (zero if not stored), binary search in the row
+	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);};
+	inline const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;};
+	const T* diagonalof(NRVec<T> &r, const bool divide=0, bool cache=false) const; //diagonal (diagonal of A^T A for non-square, as SparseMat), or divide r by it
+	const CSRMat transpose(bool conj=false) const;
+	CSRMat & operator*=(const T &a); //multiply by a scalar (the structure stays shared)
+	CSRMat & operator=(const T &a); //a times the unit matrix
+	CSRMat & operator+=(const T &a); //add a to the diagonal
+	CSRMat & operator-=(const T &a) {return *this += -a;}
+	void identity() {*this = (T)1;}
+	void axpy(const T alpha, const CSRMat &x, const bool transp=false); //this += alpha x (alpha x^T); elements cancelled to zero stay stored, see simplify()
+	CSRMat & operator+=(const CSRMat &rhs) {axpy((T)1,rhs); return *this;}
+	CSRMat & operator-=(const CSRMat &rhs) {axpy((T)-1,rhs); return *this;}
 	inline const CSRMat operator*(const T &rhs) const {return CSRMat(*this) *= rhs;}
 	inline const CSRMat operator+(const T &rhs) const {return CSRMat(*this) += rhs;}
 	inline const CSRMat operator-(const T &rhs) const {return CSRMat(*this) -= rhs;}
 	inline const CSRMat operator+(const CSRMat &rhs) const {return CSRMat(*this) += rhs;}
 	inline const CSRMat operator-(const CSRMat &rhs) const {return CSRMat(*this) -= rhs;}
-        CSRMat & operator*=(const T &a);         //multiply by a scalar
-        CSRMat & operator=(const T &a);          //assign a to diagonal
-        CSRMat & operator+=(const T &a);         //assign a to diagonal
-        CSRMat & operator-=(const T &a);         //assign a to diagonal
-	void axpy(const T alpha, const CSRMat &x, const bool transp=0); // this+= a*x
-        inline CSRMat & operator+=(const CSRMat &rhs) {axpy((T)1,rhs); return *this;};
-        inline CSRMat & operator-=(const CSRMat &rhs) {axpy((T)-1,rhs); return *this;};
-	const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal
-	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const;
-	inline const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); this->gemv((T)0,result,'n',(T)1,rhs); return result;};
-	typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
-        inline const CSRMat operator*(const CSRMat &rhs) const {CSRMat<T> r(nn,mm); r.gemm(0,*this,'n',rhs,'n',1); return r;}; 
-	void gemm(const T beta, const CSRMat &a, const char transa, const CSRMat &b, const char transb, const T alpha); //this := alpha*op( A )*op( B ) + beta*this 
-	inline SPMatindex length() {return nnz;};
-	void transposeme() const {laerror("in-place transposition not implemented for CSRMat");};
-	CSRMat transpose(bool conj=false) const {laerror("CSR transpose not implemented");}
-	void get(int fd, bool dimen, bool transp);
-        void put(int fd, bool dimen, bool transp) const;
-	int nrows() const {return nn;}
-	int ncols() const {return mm;}
-
-	/*class iterator {//not efficient, just for output to ostreams
-        private:
-                matel<T> *p;
-		matel<T> my;
-		SPMatindex mynn;
-		SPMatindex mymm;
-		
-
-        public:
-		//compiler-generated iterator & operator=(const  iterator &rhs); 
-		//compiler-generated iterator(const iterator &rhs); 
-                iterator()
-                iterator(const CSRMat &rhs)
-                iterator operator++()  
-		iterator(matel<T> *q) {p=q; col=NULL;}//just for end()
-                //compiler-generated ~iterator() {};
-                bool operator!=(const iterator &rhs) const {return p!=rhs.p;} //only for comparison with end()
-                bool operator==(const iterator &rhs) const {return p==rhs.p;} //only for comparison with end()
-                matel<T> & operator*() const {return *p;}
-                matel<T> * operator->() const {return p;}
-		bool notend() const {return (bool)p;};
-        };
-        iterator begin() const {return iterator(*this);}
-        iterator end() const {return iterator(NULL);}
-	*/
+	const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const; //Frobenius norm of A - scalar I
+	const T dot(const CSRMat &rhs) const; //\sum_{ij} conj(a_ij) b_ij as NRMat::dot
+	const T dot(const NRMat<T> &rhs) const;
+	const T trace() const;
+	void simplify(const typename LA_traits<T>::normtype sparseepsilon=LA_precision<typename LA_traits<T>::normtype>::sparse_epsilon()); //removes elements with |x| <= sparseepsilon
+	inline const NRMat<T> operator*(const NRMat<T> &rhs) const {NRMat<T> r(nn,rhs.ncols()); r.gemm((T)0,*this,'n',rhs,'n',(T)1); return r;};
+	void gemm(const T beta, const CSRMat &a, const char transa, const CSRMat &b, const char transb, const T alpha); //this = alpha op(A) op(B) + beta this, sparse result
+	inline const CSRMat operator*(const CSRMat &rhs) const {CSRMat r; r.gemm((T)0,*this,'n',rhs,'n',(T)1); return r;};
+	CSRMat & oplusequal(const CSRMat &rhs); //direct sum
+	CSRMat & oplusequal(const NRMat<T> &rhs) {return oplusequal(CSRMat(rhs));}
+	CSRMat & oplusequal(const NRSMat<T> &rhs) {return oplusequal(CSRMat(rhs));}
+	const CSRMat oplus(const CSRMat &rhs) const {return CSRMat(*this).oplusequal(rhs);}
+	const CSRMat oplus(const NRMat<T> &rhs) const {return CSRMat(*this).oplusequal(rhs);}
+	const CSRMat oplus(const NRSMat<T> &rhs) const {return CSRMat(*this).oplusequal(rhs);}
+	const CSRMat otimes(const CSRMat &rhs) const; //direct product, (A x B)_{i nB+k, j mB+l} = a_ij b_kl as NRMat::operator|
+	const CSRMat otimes(const NRMat<T> &rhs) const {return otimes(CSRMat(rhs));}
+	const CSRMat otimes(const NRSMat<T> &rhs) const {return otimes(CSRMat(rhs));}
+	//permutations: element of row r goes to row p(r) (with inverse p^{-1}); matrix indices from 0, permutation from 1; the full matrix is stored, so rows means only rows
+	void permuteme_rows(const NRPerm<int> &p, const bool inverse=false);
+	void permuteme_cols(const NRPerm<int> &p, const bool inverse=false);
+	void permuteme_both(const NRPerm<int> &p, const NRPerm<int> &q, const bool inverse=false);
+	const CSRMat permuted_rows(const NRPerm<int> &p, const bool inverse=false) const {CSRMat a(*this); a.permuteme_rows(p,inverse); return a;};
+	const CSRMat permuted_cols(const NRPerm<int> &p, const bool inverse=false) const {CSRMat a(*this); a.permuteme_cols(p,inverse); return a;};
+	const CSRMat permuted_both(const NRPerm<int> &p, const NRPerm<int> &q, const bool inverse=false) const {CSRMat a(*this); a.permuteme_both(p,q,inverse); return a;};
+	void permuterows(const NRVec<SPMatindex> &p); //as SparseMat: row r goes to row p[r], indices from 0
+	void permutecolumns(const NRVec<SPMatindex> &p);
+	void permuteindices(const NRVec<SPMatindex> &p); //rows and columns
+	void transposeme(bool conj=false) {*this=transpose(conj);}
+	void get(int fd, bool dimensions=1, bool transp=0, bool orcaformat=false);
+	void put(int fd, bool dimensions=1, bool transp=0, bool orcaformat=false) const;
 };

+template<typename T>
+inline const CSRMat<T> commutator(const CSRMat<T> &x, const CSRMat<T> &y, const bool trx=0, const bool tryy=0)
+{
+CSRMat<T> r;
+r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
+r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)-1);
+return r;
+}
+
+template<typename T>
+inline const CSRMat<T> anticommutator(const CSRMat<T> &x, const CSRMat<T> &y, const bool trx=0, const bool tryy=0)
+{
+CSRMat<T> r;
+r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
+r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)1);
+return r;
+}
+
 template <typename T>
-std::ostream & operator<<(std::ostream &s, const CSRMat<T> &x);
+inline const NRVec<T> NRVec<T>::operator*(const CSRMat<T> &mat) const
+{
+NOT_GPU(*this);
+NRVec<T> result(mat.ncols());
+result.gemv((T)0,mat,'t',(T)1,*this);
+return result;
+}

 template <class T>
-std::istream& operator>>(std::istream  &s, CSRMat<T> &x);
+std::ostream& operator<<(std::ostream &s, const CSRMat<T> &x) //the same text format as SparseMat
+		{
+		s << (int)x.nrows() << ' ' << (int)x.ncols() << '\n';
+		for(typename CSRMat<T>::iterator p=x.begin(); p.notend(); ++p)
+			s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io<T>::IOtype)p->elem << '\n';
+		s << "-1 -1\n";
+		return s;
+		}

+template <class T>
+std::istream& operator>>(std::istream &s, CSRMat<T> &x)
+		{
+		SparseMat<T> a;
+		s >> a;
+		x=CSRMat<T>(a);
+		return s;
+		}

 }//namespace
 #endif //_CSRMAT_H_
@@ -0,0 +1,2 @@
+%bibliography items of the chapter Sparse matrices (sparse.tex), cited in the comments of csrmat.cc
+\bibitem{sparse:gustavson1978} F. G. Gustavson, ACM Trans. Math. Softw. \textbf{4}, 250 (1978).
@@ -8,7 +8,7 @@ it and describe only what is specific to the class.
 \subsection{Shared data}

 An object of \code{NRVec}, \code{NRMat}, \code{NRSMat} and of all other classes holding large data
-(\code{SparseMat}, \code{SparseSMat}, \code{Tensor}, \code{fourindex}, \code{Polynomial}, the permutation
+(\code{SparseMat}, \code{SparseSMat}, \code{CSRMat}, \code{Tensor}, \code{fourindex}, \code{Polynomial}, the permutation
 classes, \code{bitvector}, \dots) consists of a small header (dimensions, a pointer to the data and a pointer
 to a reference counter) and the data block. The copy constructor and the assignment operator between objects
 of the same type copy only the header and increment the reference counter, so that both objects share the
@@ -181,7 +181,7 @@ member & double & float & used by \\
 The class templates are explicitly instantiated in the library for \code{double}, \code{std::complex<double>}
 and the integer types \code{char}, \code{short}, \code{int}, \code{long}, \code{long long} and their unsigned
 variants; \code{NRVec}, \code{NRMat}, \code{NRSMat}, \code{SparseMat}, \code{SparseSMat}, \code{Tensor}, \code{fourindex} (with \code{unsigned char} indices) and \code{Polynomial} also for \code{float} and
-\code{std::complex<float>} (single precision support, 2026, in progress for the other classes). The BLAS and LAPACK backed
+\code{std::complex<float>} (single precision support, 2026, in progress for the other classes); \code{CSRMat} only for the four floating point types. The BLAS and LAPACK backed
 functions are implemented for \code{double} and \code{std::complex<double>} (some also for mixed real--complex arguments)
 and for their single precision counterparts, which always use the single precision BLAS/LAPACK routines (\code{s}/\code{c}); calling them for other types results in a run-time
 \code{laerror} (``not implemented for this type'') or a link error. Functions implemented in the headers work
@@ -48,7 +48,7 @@ the Numerical Recipes library. They can be used with it, but are completely inde
 \header{mat.h} & \code{NRMat<T>}, \code{NRMat_from1<T>} (chapter~\ref{ch:mat}) \\
 \header{smat.h} & \code{NRSMat<T>}, \code{NRSMat_from1<T>} (chapter~\ref{ch:smat}) \\
 \header{sparsemat.h}, \header{sparsesmat.h} & sparse matrices \code{SparseMat<T>}, \code{SparseSMat<T>} (chapter~\ref{ch:sparse}) \\
-\header{csrmat.h} & compressed row storage sparse matrix (unfinished) \\
+\header{csrmat.h} & compressed row storage sparse matrix \code{CSRMat<T>} (chapter~\ref{ch:sparse}) \\
 \header{nonclass.h}, \header{matexp.h} & linear systems, eigenvalue problems, SVD, matrix functions (chapter~\ref{ch:nonclass}) \\
 \header{davidson.h}, \header{lanczos.h}, \header{conjgrad.h}, \header{gmres.h}, \header{diis.h}, \header{bisection.h} & iterative methods, DIIS convergence acceleration (chapter~\ref{ch:iterative}) \\
 \header{optimize.h} & numerical optimization (chapter~\ref{ch:optimize}) \\
@@ -111,7 +111,7 @@ its condition of use is that the vector outlives the matrix, see Pitfalls. The c
 \lafunc{explicit NRMat(const WeightPermutation<int,T> &p, const bool direction)}{Weight times permutation matrix.}
 \lafunc{explicit NRMat(const PermutationAlgebra<int,T> &p, const bool direction, const int nforce = 0)}{Sum of weighted permutation matrices of the elements of a permutation algebra; size taken from the first element unless \code{nforce>0}. (The header remarks that the permutation-matrix representation of $S_n$ is reducible to only the irreps $[n]$ and $[n-1,1]$, so e.g.\ Young projectors cannot be represented this way.) See chapter~\ref{ch:permutation}.}
 \lafunc{explicit NRMat(const SparseMat<T> &rhs); explicit NRMat(const SparseSMat<T> &rhs); explicit NRMat(const SemiSparseMat<T> &rhs)}{Dense matrix from sparse ones (chapter~\ref{ch:sparse}).}
-\lafunc{explicit NRMat(const CSRMat<T> &rhs)}{Not supported: \code{CSRMat} is unfinished, the declaration is kept for possible future work (no definition, a link error if used).}
+\lafunc{explicit NRMat(const CSRMat<T> &rhs)}{Dense matrix from a compressed row sparse one (section~\ref{sec:sparse:csr}).}
 \lafunc{~NRMat()}{Decrements the reference count, frees the data when it reaches zero.}
 \lafunc{NRMat & operator=(const NRMat &rhs)}{Shallow assignment (sec.~\ref{sec:concepts:assign}).}
 \lafunc{NRMat & operator|=(const NRMat &rhs)}{Deep copy.}
@@ -153,6 +153,7 @@ Note that scalar addition and subtraction act on the \emph{diagonal} only, i.e.\
 \lafunc{void gemm(const T &beta, const SparseMat<T> &a, const char transa, const NRMat &b, const char transb, const T &alpha)}{The same with a sparse first factor; \code{double} and \code{float} only (the complex versions call \code{laerror("not implemented yet")}), see chapter~\ref{ch:sparse}.}
 \lafunc{const NRMat operator*(const NRSMat<T> &rhs) const}{$AS$ with packed symmetric (real) or hermitian (complex) $S$ (BLAS \code{dspmv}/\code{zhpmv} per row).}
 \lafunc{SparseSMat<T> operator*(const SparseSMat<T> &rhs) const}{Dense times sparse, the result is sparse (chapter~\ref{ch:sparse}).}
+\lafunc{void gemm(const T &beta, const CSRMat<T> &a, const char transa, const NRMat &b, const char transb, const T &alpha); const NRMat operator*(const CSRMat<T> &rhs) const}{$C\leftarrow\beta C+\alpha\,\mathrm{op}(A)\,\mathrm{op}(B)$ with a compressed row sparse $A$ (all types, \code{trans} \code{'n'}, \code{'t'}, \code{'c'}), and dense times compressed row sparse with a dense result (section~\ref{sec:sparse:csr}).}
 \lafunc{const NRVec<T> operator*(const NRVec<T> &rhs) const; const NRVec<std::complex<T> > operator*(const NRVec<std::complex<T> > &rhs) const}{Matrix--vector product via \code{NRVec::gemv} (the second form: real matrix times complex vector).}
 \lafunc{void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const; void gemv(const T beta, NRVec<std::complex<T> > &r, const char trans, const T alpha, const NRVec<std::complex<T> > &x) const}{$r\leftarrow\alpha\,\mathrm{op}(A)x+\beta r$; forwards to \code{r.gemv(beta,*this,trans,alpha,x)}, see chapter~\ref{ch:vec}.}
 \lafunc{const T dot(const NRMat &rhs) const}{$\sum_{ij}A_{ij}^*B_{ij}$ (Frobenius inner product, conjugating \code{*this}).}
@@ -238,6 +239,7 @@ See section~\ref{sec:concepts:io} for the general conventions.
 Besides the constructors, \code{operator*(const SparseSMat<T>&)} and \code{gemm} with a \code{SparseMat}
 factor listed above, \code{NRMat} declares
 \begin{lareference}
+\lafunc{NRMat & operator+=(const CSRMat<T> &rhs); NRMat & operator-=(const CSRMat<T> &rhs)}{Add/subtract a compressed row sparse matrix (section~\ref{sec:sparse:csr}).}
 \lafunc{NRMat & operator+=(const SparseMat<T> &rhs); NRMat & operator-=(const SparseMat<T> &rhs)}{Add/subtract the sparse matrix (defined in \code{sparsemat.h}; a symmetric \code{SparseMat} is unfolded, without conjugation).}
 \end{lareference}

@@ -297,10 +297,10 @@ infinity for a singular matrix. It is LAPACK's estimate, a lower bound that is u
 \section{Real and imaginary parts}\label{sec:nonclass:complex}
 \begin{lareference}
 \lafunc{const LA_traits<T>::realtype realpart(const T &a); const LA_traits<T>::realtype imagpart(const T &a)}{Real or imaginary part as a real object;
-specialized for \code{NRVec}, \code{NRMat}, \code{NRSMat} and \code{SparseMat} over \code{complex<double>} and \code{complex<float>}.}
+specialized for \code{NRVec}, \code{NRMat}, \code{NRSMat}, \code{SparseMat} and \code{CSRMat} over \code{complex<double>} and \code{complex<float>}.}
 \lafunc{const LA_traits<T>::complextype realmatrix(const T &a); const LA_traits<T>::complextype imagmatrix(const T &a)}{Complex object $A$ resp.\ $\mathrm{i}A$ from a real one;
-for \code{NRMat} and \code{SparseSMat} of \code{double} and \code{float}.}
-\lafunc{const LA_traits<T>::complextype complexmatrix(const T &re, const T &im)}{Complex object $\mathrm{re}+\mathrm{i}\,\mathrm{im}$; for \code{NRMat} and \code{SparseSMat} of \code{double} and \code{float}.}
+for \code{NRMat}, \code{SparseSMat} and \code{CSRMat} of \code{double} and \code{float}.}
+\lafunc{const LA_traits<T>::complextype complexmatrix(const T &re, const T &im)}{Complex object $\mathrm{re}+\mathrm{i}\,\mathrm{im}$; for \code{NRMat}, \code{SparseSMat} and \code{CSRMat} of \code{double} and \code{float}.}
 \end{lareference}
 For the same purpose the classes provide \code{complexify()} (chapters~\ref{ch:vec}, \ref{ch:mat}).
 The real and imaginary part functions are CPU only (sec.~\ref{sec:concepts:gpu}).
@@ -107,7 +107,7 @@ are CPU only (\code{NOT_GPU}).
 \lafunc{NRSMat(const typename LA_traits_complex<T>::NRSMat_Noncomplex_type &rhs, bool imagpart = false)}{Complex matrix from a real one: \code{rhs} becomes the real part (\code{imagpart=false}) or the imaginary part (\code{imagpart=true}), the other part is zero. Defined only for \code{T=std::complex<double>} and \code{std::complex<float>}.}
 \lafunc{explicit NRSMat(const NRMat<T> &rhs)}{Symmetric part $(A+A\Tr)/2$ of a square dense matrix (deep copy; nonsquare input is checked only in DEBUG).}
 \lafunc{explicit NRSMat(const NRVec<T> &rhs, const int n)}{Interprets a vector of length $n(n+1)/2$ as packed matrix of order $n$. \emph{Shares} the data with \code{rhs} (no copy); length checked only in DEBUG.}
-\lafunc{explicit NRSMat(const SparseMat<T> &rhs); explicit NRSMat(const SparseSMat<T> &rhs)}{Dense from sparse, see \ref{ch:sparse}.}
+\lafunc{explicit NRSMat(const SparseMat<T> &rhs); explicit NRSMat(const SparseSMat<T> &rhs); explicit NRSMat(const CSRMat<T> &rhs)}{Dense from sparse, see \ref{ch:sparse}; from \code{SparseSMat} and \code{CSRMat} the lower triangle is taken.}
 \lafunc{~NRSMat()}{Destructor, frees the data when the last reference disappears.}
 \lafunc{NRSMat & operator=(const NRSMat &rhs)}{Shallow copy, see section~\ref{sec:concepts:assign}.}
 \lafunc{NRSMat & operator|=(const NRSMat &rhs)}{Deep copy.}
@@ -1,8 +1,8 @@
 \chapter{Sparse matrices}\label{ch:sparse}
-\laheader{sparsemat.h, sparsemat.cc, sparsesmat.h, sparsesmat.cc}
+\laheader{sparsemat.h, sparsemat.cc, sparsesmat.h, sparsesmat.cc, csrmat.h, csrmat.cc}

 \section{Overview}\label{sec:sparse:overview}
-\LA{} offers two independent sparse matrix classes with different storage schemes, plus a small helper class:
+\LA{} offers three independent sparse matrix classes with different storage schemes, plus a small helper class:
 \begin{itemize}
 \item \code{SparseMat<T>} (section~\ref{sec:sparse:sparsemat}): an \emph{unordered linked list} of
  (row, column, value) triplets. Appending elements is cheap, duplicates are allowed and summed,
@@ -15,10 +15,12 @@
  exponentiating symmetric matrices (e.g.\ the unitary complex symmetric $\exp(\mathrm{i}H)$ with real symmetric $H$),
  but despite its name it is a general (also nonsymmetric and rectangular) sparse matrix; symmetry
  (hermiticity for complex types) is only a convention of some of its operations.
+\item \code{CSRMat<T>} (section~\ref{sec:sparse:csr}): the compressed sparse row format (row pointers, column
+  indices, values), read-efficient for products with vectors and matrices, e.g.\ in iterative solvers; built by
+  conversion from the other classes.
 \item \code{SemiSparseMat<T>} (section~\ref{sec:sparse:semisparse}): a dense diagonal plus a
  \code{SparseMat} off-diagonal part, just enough to be used in iterative eigensolvers.
 \end{itemize}
-\code{CSRMat<T>} (section~\ref{sec:sparse:csr}) is unfinished.

 Common definitions in \code{sparsemat.h}:
 \begin{lareference}
@@ -88,6 +90,7 @@ Most modifying members call \code{copyonwrite()}: \code{simplify}, \code{operato
 \lafunc{SparseMat(const SparseMat &rhs)}{Shallow copy; an empty matrix gets its own reference count.}
 \lafunc{explicit SparseMat(const NRMat<T> &rhs)}{From a dense matrix, elements with $|a_{ij}|>$ $\epsilon_{\rm sp}$.}
 \lafunc{explicit SparseMat(const NRSMat<T> &rhs)}{From a packed symmetric matrix; symmetric flag set, lower triangle stored.}
+\lafunc{explicit SparseMat(const CSRMat<T> &rhs)}{From a \code{CSRMat} (section~\ref{sec:sparse:csr}): all stored elements, not symmetric.}
 \lafunc{~SparseMat()}{Destructor.}
 \lafunc{SparseMat & operator=(const SparseMat &rhs)}{Shallow copy (also copies the symmetric flag).}
 \lafunc{SparseMat & operator=(const T &a)}{Scalar matrix $a\cdot1$; the matrix must be square and dimensioned, old contents are discarded.}
@@ -248,7 +251,7 @@ copies).
 \lafunc{explicit SparseSMat(const NRSMat<T> &rhs)}{From a packed symmetric (complex: hermitian) matrix (both triangles inserted, the upper one conjugated, $|a_{ij}|>$ $\epsilon_{\rm sp}$).}
 \lafunc{explicit SparseSMat(const NRMat<T> &rhs)}{From a dense $n\times m$ matrix, nonzero elements copied as they are (no symmetrization).}
 \lafunc{explicit SparseSMat(const SparseMat<T> &rhs)}{From a \code{SparseMat} (same dimensions); duplicates are summed, a symmetric \code{SparseMat} is unfolded without conjugation (complex symmetric, as \code{SparseMat} defines it).}
-\lafunc{explicit SparseSMat(const CSRMat<T> &rhs)}{Declared but not implemented (\code{CSRMat} is not supported).}
+\lafunc{explicit SparseSMat(const CSRMat<T> &rhs)}{From a \code{CSRMat} (section~\ref{sec:sparse:csr}): all stored elements, without mirroring.}
 \lafunc{~SparseSMat()}{Destructor.}
 \lafunc{SparseSMat & operator=(const SparseSMat &rhs)}{Shallow copy.}
 \lafunc{SparseSMat & operator=(const T &a)}{Scalar matrix $a\cdot1$ (diagonal elements are created even for $a=0$).}
@@ -310,11 +313,100 @@ In \ref{ch:nonclass}: \code{complexmatrix()}, \code{realmatrix()} and \code{imag
 \code{SparseSMat}.

 \section{CSRMat}\label{sec:sparse:csr}
-\laheader{csrmat.h}
-\code{CSRMat<T>} is meant to hold a matrix in compressed sparse row format (row pointers, column indices,
-values, optionally on the GPU via cuSPARSE) for fast products. The class is unfinished: only the declarations
-exist, \code{csrmat.cc} contains no implementation and no instantiation. It is \emph{not supported} and must not be used;
-the declarations (also \code{NRMat(const CSRMat&)} and \code{SparseSMat(const CSRMat&)}) are kept for possible future work.
+\laheader{csrmat.h, csrmat.cc}
+\code{CSRMat<T>} stores a matrix in the compressed sparse row format: the elements of row $i$ are \code{v[k]} in
+column \code{col[k]} for $k=$ \code{rowptr[i]}$,\ldots,$\code{rowptr[i+1]}$-1$, with the columns strictly ascending within
+each row (an invariant of the class). The format is read-efficient: $Ax$ is a row-wise dot product without any
+indirection in the result, which makes it the fastest of the sparse classes in iterative solvers (\ref{ch:iterative})
+and matrix functions with a fixed sparsity pattern. It cannot be built element by element; a matrix is assembled as a
+\code{SparseMat} (\code{add}) or a \code{SparseSMat} (or as three arrays) and converted, in $O(n_{\rm nz}+n+m)$.
+
+The full matrix is always stored, also for symmetric and hermitian matrices. The flag \code{symmetric} only
+\emph{declares} the matrix hermitian, $A=A^\dagger$ (real symmetric for real types); it is not checked and never changes the
+arithmetic, but it is what \code{issymmetric()} returns to the eigensolvers (\code{davidson} uses the hermitian algorithm,
+\code{lanczos} requires it). This is the convention of \code{NRSMat} and \code{SparseSMat}, not the complex symmetric
+one of \code{SparseMat}.
+
+The three arrays are \code{NRVec} members (\code{rowptr} of length $n+1$, \code{col} and \code{v} of length $n_{\rm nz}$, the
+number of stored elements); each of them shares its data independently (section~\ref{sec:concepts:refcount}). Operations which
+change only the values keep the index arrays shared, e.g.\ \code{B=A*2.}, \code{A+=a} with all diagonal elements stored, or
+\code{realmatrix(A)}; \code{copyonwrite()} copies all three. Conversions store only elements with $|a_{ij}|>\epsilon_{\rm sp}$,
+whereas arithmetic keeps elements which cancel to zero, until \code{simplify()} is called. $n_{\rm nz}$ is limited to $2^{31}-1$
+(the \code{NRVec} length). The class is instantiated for \code{double}, \code{float}, \code{std::complex<double>} and
+\code{std::complex<float>}; all operations are CPU only (\code{laerror} for vectors and matrices on the GPU), the device-dependent
+format is the one of cuSPARSE, so a GPU version can be added later.
+
+\subsection{Construction, conversion, access}
+\begin{lareference}
+\lafunc{CSRMat()}{Empty $0\times0$ matrix.}
+\lafunc{explicit CSRMat(const SPMatindex n, const SPMatindex m); explicit CSRMat(const SPMatindex n)}{Zero $n\times m$ resp.\ $n\times n$ matrix.}
+\lafunc{CSRMat(const SPMatindex n, const SPMatindex m, const NRVec<SPMatindex> &rowptr, const NRVec<SPMatindex> &col, const NRVec<T> &v)}{From the three arrays, which are shared (not copied). In DEBUG their consistency, the column range and the ascending columns are checked.}
+\lafunc{explicit CSRMat(const SparseMat<T> &rhs)}{Duplicates are summed, elements with $|a_{ij}|\le\epsilon_{\rm sp}$ dropped. A symmetric \code{SparseMat} is unfolded without conjugation (complex symmetric, as \code{SparseMat} defines it); the flag is set for the real types only.}
+\lafunc{explicit CSRMat(const SparseSMat<T> &rhs); explicit CSRMat(const NRMat<T> &rhs)}{Elements with $|a_{ij}|>\epsilon_{\rm sp}$; the flag is not set (call \code{setsymmetric()} if appropriate).}
+\lafunc{explicit CSRMat(const NRSMat<T> &rhs)}{Hermitian as \code{NRSMat}: the upper triangle is the conjugate of the stored lower one; the flag is set.}
+\lafunc{void copyonwrite(bool detachonly=false, bool deep=true); void resize(const SPMatindex n, const SPMatindex m); void dealloc(); void clear(); void swap(CSRMat &rhs)}{As for the other classes; \code{resize} gives a zero matrix and clears the flag, \code{clear} keeps it.}
+\lafunc{SPMatindex nrows() const; SPMatindex ncols() const; int length() const}{Dimensions and $n_{\rm nz}$.}
+\lafunc{bool issymmetric() const; void setsymmetric(); void setunsymmetric()}{The hermitian flag; \code{setsymmetric} only declares (square matrices only).}
+\lafunc{const NRVec<SPMatindex> & getrowptr() const; const NRVec<SPMatindex> & getcol() const; const NRVec<T> & getval() const}{Read access to the arrays.}
+\lafunc{const T operator()(const SPMatindex i, const SPMatindex j) const}{Element $a_{ij}$ by binary search in row $i$, zero if not stored. There is no write access to single elements.}
+\lafunc{iterator begin() const; iterator end() const}{Read-only iterator over the stored elements in row-major order: \code{p->row}, \code{p->col}, \code{p->elem} (a \code{matel<T>}), \code{p.notend()}.}
+\end{lareference}
+Other classes from \code{CSRMat}: \code{explicit NRMat<T>(const CSRMat<T>&)} (full matrix), \code{explicit NRSMat<T>(const CSRMat<T>&)}
+(lower triangle, square only, as from \code{SparseSMat}), \code{explicit SparseMat<T>(const CSRMat<T>&)} (not symmetric),
+\code{explicit SparseSMat<T>(const CSRMat<T>&)} and \code{explicit NRVec<T>(const CSRMat<T>&)} (for a $1\times m$ or $n\times 1$ matrix).
+
+\subsection{Arithmetic}
+\begin{lareference}
+\lafunc{CSRMat & operator*=(const T &a); const CSRMat operator*(const T &a) const}{Scaling of the values (BLAS); the index arrays stay shared, $a=0$ gives the zero matrix.}
+\lafunc{CSRMat & operator=(const T &a); void identity()}{$a\cdot1$ resp.\ the unit matrix (square only).}
+\lafunc{CSRMat & operator+=(const T &a); CSRMat & operator-=(const T &a); const CSRMat operator+(const T &a) const; const CSRMat operator-(const T &a) const}{$A\pm a\cdot1$ (square only): in place if all diagonal elements are stored, otherwise by a merge.}
+\lafunc{void axpy(const T alpha, const CSRMat &x, const bool transp=false)}{$A\leftarrow A+\alpha x$ (or $\alpha x\Tr$): merge of the sorted rows in two passes; for the same structure (shared or equal index arrays) a BLAS \code{axpy} of the values. Elements which cancel stay stored.}
+\lafunc{CSRMat & operator+=(const CSRMat &rhs); CSRMat & operator-=(const CSRMat &rhs); const CSRMat operator+(const CSRMat &rhs) const; const CSRMat operator-(const CSRMat &rhs) const}{Via \code{axpy}.}
+\lafunc{const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const}{Frobenius norm of $A-\mathrm{scalar}\cdot1$.}
+\lafunc{const T dot(const CSRMat &rhs) const; const T dot(const NRMat<T> &rhs) const}{$\sum_{ij}\bar a_{ij}b_{ij}$, conjugating the left operand as \code{NRMat::dot} (unlike \code{SparseMat::dot}).}
+\lafunc{const T trace() const}{Sum of the diagonal elements.}
+\lafunc{void simplify(const typename LA_traits<T>::normtype sparseepsilon=LA_precision<normtype>::sparse_epsilon())}{Removes the elements with $|a_{ij}|\le$ \code{sparseepsilon}.}
+\end{lareference}
+The result of an addition is declared hermitian if both operands are (a matrix without stored elements counts as hermitian) and the
+scalar factor is real; multiplying by or adding a non-real scalar clears the flag. In other classes:
+\code{NRMat<T>::operator+=(const CSRMat<T>&)} and \code{operator-=}.
+
+\subsection{Products}
+\begin{lareference}
+\lafunc{void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const}{$r\leftarrow\beta r+\alpha\,\mathrm{op}(A)x$, \code{trans} \code{'n'}, \code{'t'} or \code{'c'} (conjugate transposition, the same as \code{'t'} for real types); forwards to \code{NRVec<T>::gemv}. The size of \code{x} and (in DEBUG) of the result are checked, $\beta=0$ clears the result.}
+\lafunc{inline const NRVec<T> operator*(const NRVec<T> &rhs) const}{$Ax$.}
+\lafunc{const T* diagonalof(NRVec<T> &r, const bool divide=0, bool cache=false) const}{Diagonal (binary search in each row), or $r$ divided by it elementwise (zero diagonal elements skipped). For a non-square matrix the diagonal of $A\Tr A$ ($\sum_i a_{ij}^2$, no conjugation, as \code{SparseMat}), used as the preconditioner of \code{conjgrad}. \code{cache} is ignored.}
+\lafunc{inline const NRMat<T> operator*(const NRMat<T> &rhs) const}{$AB$ with a dense result, via \code{NRMat<T>::gemm} below.}
+\lafunc{void gemm(const T beta, const CSRMat &a, const char transa, const CSRMat &b, const char transb, const T alpha)}{$C\leftarrow\beta C+\alpha\,\mathrm{op}(A)\,\mathrm{op}(B)$ with a sparse result (\code{trans} as for \code{gemv}), by the row-wise algorithm of Gustavson with a dense accumulator \cite{sparse:gustavson1978} in two passes; the rows of the result are sorted, elements which cancel stay stored, the result is not declared hermitian. With $\beta=0$ the dimensions of \code{*this} do not matter.}
+\lafunc{inline const CSRMat operator*(const CSRMat &rhs) const}{$AB$.}
+\end{lareference}
+In other classes: \code{NRVec<T>::gemv(const T beta, const CSRMat<T> &a, const char trans, const T alpha, const NRVec &x)},
+\code{NRVec<T>::operator*(const CSRMat<T>&)} ($x\Tr A$),
+\code{NRMat<T>::gemm(const T &beta, const CSRMat<T> &a, const char transa, const NRMat &b, const char transb, const T &alpha)}
+(all combinations of \code{'n'}, \code{'t'}, \code{'c'}) and \code{NRMat<T>::operator*(const CSRMat<T>&)} (dense times sparse, dense result).
+The sparse kernels are plain loops (there is no sparse BLAS in the libraries used).
+
+\subsection{Transposition, permutations, direct sums and products}
+\begin{lareference}
+\lafunc{const CSRMat transpose(bool conj=false) const; void transposeme(bool conj=false)}{$A\Tr$ or $A\Hc$ by a counting sort, $O(n_{\rm nz}+n+m)$; the flag is kept.}
+\lafunc{void permuteme_rows(const NRPerm<int> &p, const bool inverse=false); void permuteme_cols(const NRPerm<int> &p, const bool inverse=false); void permuteme_both(const NRPerm<int> &p, const NRPerm<int> &q, const bool inverse=false)}{The convention of \code{SparseMat}: element $(i,j)$ moves to row $p_{i+1}-1$ resp.\ column $q_{j+1}-1$ (matrix indices from 0, permutations from 1, optionally inverted). Since the full matrix is stored, \code{permuteme_rows} permutes only the rows also of a hermitian matrix (unlike a symmetric \code{SparseMat}); the flag is cleared, except by \code{permuteme_both} with $p=q$.}
+\lafunc{const CSRMat permuted_rows(...) const; const CSRMat permuted_cols(...) const; const CSRMat permuted_both(...) const}{Out-of-place versions.}
+\lafunc{void permuterows(const NRVec<SPMatindex> &p); void permutecolumns(const NRVec<SPMatindex> &p); void permuteindices(const NRVec<SPMatindex> &p)}{As for \code{SparseMat}: $i\to p_i$ (indices from 0), \code{permuteindices} for rows and columns.}
+\lafunc{CSRMat & oplusequal(const CSRMat &rhs); const CSRMat oplus(const CSRMat &rhs) const}{Direct sum $A\oplus B$; also with \code{NRMat} and \code{NRSMat} (converted as above, i.e.\ \code{NRSMat} as hermitian).}
+\lafunc{const CSRMat otimes(const CSRMat &rhs) const}{Direct (Kronecker) product, $(A\otimes B)_{i n_B+k,\,j m_B+l}=a_{ij}b_{kl}$ as \code{NRMat::operator|}; also with \code{NRMat} and \code{NRSMat}.}
+\end{lareference}
+The results of \code{oplus} and \code{otimes} are declared hermitian if both operands are.
+
+\subsection{Input/output and free functions}
+\begin{lareference}
+\lafunc{void put(int fd, bool dimensions=1, bool transp=0, bool orcaformat=false) const; void get(int fd, bool dimensions=1, bool transp=0, bool orcaformat=false)}{Binary I/O: if \code{dimensions}, $n$ and $m$ (\code{SPMatindex}); the flag (\code{int}); the three arrays by \code{NRVec::put} with their lengths. \code{transp} writes or reads the transposed matrix; \code{orcaformat=true} raises \code{laerror}.}
+\lafunc{std::ostream& operator<<(std::ostream &s, const CSRMat<T> &x); std::istream& operator>>(std::istream &s, CSRMat<T> &x)}{Text I/O in the format of \code{SparseMat} (read through a \code{SparseMat}, so duplicates are summed); the flag is not written.}
+\lafunc{inline const CSRMat<T> commutator(const CSRMat<T> &x, const CSRMat<T> &y, const bool trx=0, const bool tryy=0); inline const CSRMat<T> anticommutator(...)}{$xy\mp yx$ via two \code{gemm} calls.}
+\end{lareference}
+In \ref{ch:nonclass}: \code{realpart()}, \code{imagpart()} of a complex \code{CSRMat} and \code{realmatrix()}, \code{imagmatrix()}, \code{complexmatrix()}
+of real ones; only the values are converted (BLAS), the index arrays stay shared unless elements which became zero are removed.
+The real part of a hermitian matrix and \code{realmatrix} of a declared symmetric one keep the flag. The templates of \code{matexp.h}
+(\code{exp()}, \code{sincos()}, \code{exptimes()}) and the iterative solvers of \ref{ch:iterative} work with \code{CSRMat}.

 \section{Examples}\label{sec:sparse:examples}
 \begin{lstlisting}
@@ -353,6 +445,15 @@ h.add(2,2,1.);
 SparseSMat<double> h2 = h*h;             //result stays symmetric
 NRVec<double> w = h*x;
 std::cout << y << z << w << NRMat<double>(h2) << sd;
+
+//compressed row format: assembled as SparseMat, then converted
+SparseMat<double> hb(4,4);
+hb.add(0,0,1.); hb.add(1,0,0.2); hb.add(0,1,0.2); hb.add(3,3,2.);
+CSRMat<double> hc(hb);
+hc.setsymmetric();                       //declare it symmetric, e.g. for davidson
+NRVec<double> u = hc*x;
+CSRMat<std::complex<double> > ev = exp(imagmatrix(hc)); //exp(iH), sparse
+CSRMat<double> hc2 = hc*2.;              //shares rowptr and col with hc
 }
 \end{lstlisting}

@@ -364,4 +465,7 @@ std::cout << y << z << w << NRMat<double>(h2) << sd;
 \item In a symmetric \code{SparseMat} never store both $a_{ij}$ and $a_{ji}$; \code{setsymmetric()} does not check this.
 \item \code{SparseSMat::gemm} computes $A\Tr B$ (or $A\Hc B$); use \code{operator*} for the general product.
 \item Complex \code{SparseMat} matrices are complex symmetric, never hermitian, in all mirroring operations, whereas \code{SparseSMat} (\code{add} with \code{both=true}, conversion from \code{NRSMat}, \code{otimes_sparse}) is hermitian.
+\item \code{CSRMat} cannot be filled element by element; assemble a \code{SparseMat} or \code{SparseSMat} and convert it.
+\item The flag of \code{CSRMat} only declares hermiticity. Conversions from \code{NRMat}, \code{SparseSMat} and arrays (and from a complex
+  symmetric \code{SparseMat}) do not set it; call \code{setsymmetric()} for a hermitian matrix before \code{davidson} or \code{lanczos}.
 \end{itemize}
@@ -66,6 +66,7 @@ one: \code{rhs} becomes the real part (or the imaginary part if \code{imagpart})
 $n(n+1)/2$ elements. \emph{Shares} the data with \code{S} (no copy).}
 \lafunc{explicit NRVec(const NRMat<T> &rhs)}{View of the row-major storage of an $n\times m$ matrix as a vector of $nm$ elements.
 \emph{Shares} the data with \code{rhs}; if the library is compiled with \code{MATPTR}, the data are copied instead.}
+\lafunc{explicit NRVec(const CSRMat<T> &rhs)}{The same from a compressed row sparse matrix (section~\ref{sec:sparse:csr}).}
 \lafunc{explicit NRVec(const SparseMat<T> &rhs)}{Dense vector from a sparse matrix having one row or one column (otherwise
 \code{laerror}); duplicate entries are summed. Defined in \texttt{sparsemat.cc} for \code{double} and \code{std::complex<double>}.}
 \lafunc{~NRVec()}{Decrements the reference count, frees the data when it reaches zero.}
@@ -162,6 +163,7 @@ with $y=$ \code{*this}; \code{trans} is \code{'n'} ($A$), \code{'t'} ($A\Tr$) or
 \code{*this} must already have the correct size (checked in DEBUG builds) and defined contents if $\beta\neq0$.}
 \lafunc{void gemv(const T beta, const NRSMat<T> &a, const char trans, const T alpha, const NRVec &x)}{Same for a packed
 symmetric (real) or Hermitian (complex) matrix (BLAS spmv/hpmv); \code{trans}: for complex \code{'t'} gives $a\Tr x = a^* x$, \code{'n'} and \code{'c'} give $a x$; ignored for real matrices.}
+\lafunc{void gemv(const T beta, const CSRMat<T> &a, const char trans, const T alpha, const NRVec &x)}{The same for a compressed row sparse matrix, \code{trans} \code{'n'}, \code{'t'} or \code{'c'} (section~\ref{sec:sparse:csr}).}
 \lafunc{void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x, const bool treat_as_symmetric = false)}{Same
 for a sparse matrix (defined in \texttt{sparsemat.cc}, see Chapter~\ref{ch:sparse}); with \code{treat_as_symmetric} only one
 triangle is assumed to be stored.}
@@ -172,7 +174,7 @@ real and imaginary parts are processed separately.}
 \lafunc{const NRVec operator*(const NRMat<T> &mat) const}{Row vector times matrix, $x\Tr A$ (gemv with \code{'t'}; no conjugation).
 For $A x$ see \code{NRMat::operator*} in Chapter~\ref{ch:mat}.}
 \lafunc{const NRVec operator*(const NRSMat<T> &mat) const}{Row vector times packed symmetric/hermitian matrix, $x\Tr S$ (for hermitian $S$ computed as $S\Tr x = S^* x$ by \code{gemv} with \code{'t'}).}
-\lafunc{const NRVec operator*(const SparseMat<T> &mat) const}{Row vector times sparse matrix; CPU only.}
+\lafunc{const NRVec operator*(const SparseMat<T> &mat) const; const NRVec operator*(const CSRMat<T> &mat) const}{Row vector times sparse matrix; CPU only.}
 \lafunc{const NRMat<T> otimes(const NRVec<T> &rhs, const bool conjugate = false, const T &scale = 1) const}{Outer product
 $A_{ij}=\mathrm{scale}\cdot x_i\, y_j$, with $y_j^*$ if \code{conjugate} (complex only); BLAS ger/geru/gerc.}
 \lafunc{const NRMat<T> operator|(const NRVec<T> &rhs) const}{Same as \code{otimes(rhs,true)}, i.e.\ $x y\Hc$.}
@@ -424,9 +424,13 @@ public:
 	NRMat & operator+=(const SparseMat<T> &rhs);
 	//! subtract given sparse matrix
        NRMat & operator-=(const SparseMat<T> &rhs);
+	NRMat & operator+=(const CSRMat<T> &rhs);
+	NRMat & operator-=(const CSRMat<T> &rhs);

 	//! perform the \b gemm operation
 	void gemm(const T &beta, const SparseMat<T> &a, const char transa, const NRMat &b, const char transb, const T &alpha);
+	void gemm(const T &beta, const CSRMat<T> &a, const char transa, const NRMat &b, const char transb, const T &alpha);
+	const NRMat operator*(const CSRMat<T> &rhs) const; //dense times compressed row sparse, dense result

        inline void simplify() {};
 	bool issymmetric() const { return 0; };
@@ -235,6 +235,7 @@ public:
 //members concerning sparse matrix
 	explicit NRSMat(const SparseMat<T> &rhs);               // dense from sparse
 	explicit NRSMat(const SparseSMat<T> &rhs);               // dense from sparse
+	explicit NRSMat(const CSRMat<T> &rhs);               // dense from sparse (lower triangle)
 	inline void simplify() {}; //just for compatibility with sparse ones
 	bool issymmetric() const {return 1;}
 	void swap(NRSMat &rhs) //more efficient swap than via tmp and constructors and operator=
@@ -80,6 +80,7 @@ public:
 	inline int getcount() const {return count?*count:0;}
 	explicit SparseMat(const NRMat<T> &rhs); //construct from a dense one
 	explicit SparseMat(const NRSMat<T> &rhs); //construct from a dense symmetric one
+	explicit SparseMat(const CSRMat<T> &rhs); //full storage, not symmetric
 	SparseMat & operator=(const SparseMat &rhs);
 	SparseMat & operator=(const T &a);          //assign a to diagonal
 	void identity() {*this = (T)1;};
@@ -1880,7 +1880,6 @@ c.moveto(gpu1);
 cout << "Error = "<<(c-cgpu).norm()<<endl;
 }

-/*
 if(0)
 {
 CSRMat<double> h0;
@@ -1903,7 +1902,6 @@ NRMat<complex<double> >r2 = exp(h2);
 cout <<"error = "<<(r2-NRMat<complex<double> >(r)).norm()<<endl;
 }
 }
-*/


 if(0)
@@ -337,6 +337,7 @@ public:
 	void gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec &x);
 	void gemv(const T beta, const NRSMat<T> &a, const char trans /**< just for compatibility reasons */, const T alpha, const NRVec &x);
 	void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x,const bool treat_as_symmetric = false);
+	void gemv(const T beta, const CSRMat<T> &a, const char trans, const T alpha, const NRVec &x);

 	void gemv(	const typename LA_traits_complex<T>::Component_type beta, 
 			const typename LA_traits_complex<T>::NRMat_Noncomplex_type &a, 
@@ -376,6 +377,9 @@ public:
 		return result;
 	};

+	//! multiply given compressed row sparse matrix with this vector from left (defined in csrmat.h)
+	inline const NRVec operator*(const CSRMat<T> &mat) const;
+	
 	//! compute the outer product of two vectors 
 	const NRMat<T> otimes(const NRVec<T> &rhs, const bool conjugate = false, const T &scale = 1) const;

@@ -486,6 +490,7 @@ public:

 	//! constructor creating vector from sparse matrix
 	explicit NRVec(const SparseMat<T> &rhs);
+	explicit NRVec(const CSRMat<T> &rhs);
 	
 	//! routine for compatibility with sparse types
 	inline void simplify() {};