*** empty log message ***

mat.cc

@@ -14,43 +14,7 @@ template NRMat<char>;
  *  Templates first, specializations for BLAS next
  */
 
-// dtor
-template <typename T>
-NRMat<T>::~NRMat()
-{
-	if (!count) return;
-	if (--(*count) <= 0) {
-		if (v) {
-#ifdef MATPTR
-			delete[] (v[0]);
-#endif
-			delete[] v;
-		}
-		delete count;
-	}
-}
 
-// assign NRMat = NRMat
-template <typename T>
-NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
-{
-	if (this == &rhs) return *this;
-	if (count) {
-		if (--(*count) ==0 ) {
-#ifdef MATPTR
-			delete[] (v[0]);
-#endif
-			delete[] v;
-			delete count;
-		}
-		v = rhs.v;
-		nn = rhs.nn;
-		mm = rhs.mm;
-		count = rhs.count;
-		if (count) (*count)--;
-	}
-	return *this;
-}
 
 // Assign diagonal
 template <typename T>
@@ -68,54 +32,7 @@ NRMat<T> & NRMat<T>::operator=(const T &a)
 	 return *this;
 }
 
-// Explicit deep copy of NRmat
-template <typename T>
-NRMat<T> & NRMat<T>::operator|=(const NRMat<T> &rhs)
-{
-	if (this == &rhs) return *this;
-#ifdef DEBUG
-	if (!rhs.v) laerror("unallocated rhs in Mat operator |=");
-#endif
-	if (count)
-		if (*count > 1) {
-			--(*count);
-			nn = 0;
-			mm = 0;
-			count = 0;
-			v = 0;
-		}
-	if (nn != rhs.nn || mm != rhs.mm) {
-		if (v) {
-#ifdef MATPTR
-			delete[] (v[0]);
-#endif
-			delete[] (v);
-			v = 0;
-		}
-		nn = rhs.nn;
-		mm = rhs.mm;
-	}
-	if (!v) {
-#ifdef MATPTR
-		v = new T*[nn];
-		v[0] = new T[mm*nn];
-#else
-		v = new T[mm*nn];
-#endif
-	}
 
-#ifdef MATPTR
-	for (int i=1; i< nn; i++) v[i] = v[i-1] + mm;
-	memcpy(v[0], rhs.v[0], nn*mm*sizeof(T));
-#else
-	memcpy(v, rhs.v, nn*mm*sizeof(T));
-#endif
-
-	if (!count) count = new int;
-	*count = 1;
-
-	return *this;
-}
 
 // M += a
 template <typename T>
@@ -215,76 +132,7 @@ const NRVec<T> NRMat<T>::rsum() const
 	return result;
 }
 
-// make detach Mat and make it's own deep copy
-template <typename T>
-void NRMat<T>::copyonwrite()
-{
-#ifdef DEBUG
-	if (!count) laerror("Mat::copyonwrite of undefined matrix");
-#endif
-	if (*count > 1) {
-		(*count)--;
-		count = new int;
-		*count = 1;
-#ifdef MATPTR
-		T **newv = new T*[nn];
-		newv[0] = new T[mm*nn];
-		memcpy(newv[0], v[0], mm*nn*sizeof(T));
-		v = newv;
-		for (int i=1; i< nn; i++) v[i] = v[i-1] + mm;
-#else
-		T *newv = new T[mm*nn];
-		memcpy(newv, v, mm*nn*sizeof(T));
-		v = newv;
-#endif
-	}
-}
 
-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()");
-#endif
-	if (count)
-		if (*count > 1) {
-			(*count)--;
-			count = 0;
-			v  = 0;
-			nn = 0;
-			mm = 0;
-		}
-	if (!count) {
-		count = new int;
-		*count = 1;
-		nn = n;
-		mm = m;
-#ifdef MATPTR
-		v = new T*[nn];
-		v[0] = new T[m*n];
-		for (int i=1; i< n; i++) v[i] = v[i-1] + m;
-#else
-		v = new T[m*n];
-#endif
-		return;
-	}
-	// At this point *count = 1, check if resize is necessary
-	if (n!=nn || m!=mm) {
-		nn = n;
-		mm = m;
-#ifdef MATPTR
-		delete[] (v[0]);
-#endif
-		delete[] v;
-#ifdef MATPTR
-		v = new T*[nn];
-		v[0] = new T[m*n];
-		for (int i=1; i< n; i++) v[i] = v[i-1] + m;
-#else
-		v = new T[m*n];
-#endif
-	}
-}
 
 // transpose Mat
 template <typename T>

mat.h

@@ -327,6 +327,173 @@ inline const complex<double>  NRMat< complex<double> >::amax() const
 }
 
 
+//basi stuff to be available for any type ... must be in .h
+// dtor
+template <typename T>
+NRMat<T>::~NRMat()
+{
+        if (!count) return;
+        if (--(*count) <= 0) {
+                if (v) {
+#ifdef MATPTR
+                        delete[] (v[0]);
+#endif
+                        delete[] v;
+                }
+                delete count;
+        }
+}
+
+// assign NRMat = NRMat
+template <typename T>
+NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
+{
+        if (this == &rhs) return *this;
+        if (count) {
+                if (--(*count) ==0 ) {
+#ifdef MATPTR
+                        delete[] (v[0]);
+#endif
+                        delete[] v;
+                        delete count;
+                }
+                v = rhs.v;
+                nn = rhs.nn;
+                mm = rhs.mm;
+                count = rhs.count;
+                if (count) (*count)--;
+        }
+        return *this;
+}
+
+
+// Explicit deep copy of NRmat
+template <typename T>
+NRMat<T> & NRMat<T>::operator|=(const NRMat<T> &rhs)
+{
+        if (this == &rhs) return *this;
+#ifdef DEBUG
+        if (!rhs.v) laerror("unallocated rhs in Mat operator |=");
+#endif
+        if (count)
+                if (*count > 1) {
+                        --(*count);
+                        nn = 0;
+                        mm = 0;
+                        count = 0;
+                        v = 0;
+                }
+        if (nn != rhs.nn || mm != rhs.mm) {
+                if (v) {
+#ifdef MATPTR
+                        delete[] (v[0]);
+#endif
+                        delete[] (v);
+                        v = 0;
+                }
+                nn = rhs.nn;
+                mm = rhs.mm;
+        }
+        if (!v) {
+#ifdef MATPTR
+                v = new T*[nn];
+                v[0] = new T[mm*nn];
+#else
+                v = new T[mm*nn];
+#endif
+        }
+#ifdef MATPTR
+        for (int i=1; i< nn; i++) v[i] = v[i-1] + mm;
+        memcpy(v[0], rhs.v[0], nn*mm*sizeof(T));
+#else
+        memcpy(v, rhs.v, nn*mm*sizeof(T));
+#endif
+
+        if (!count) count = new int;
+        *count = 1;
+
+        return *this;
+}
+
+// make detach Mat and make it's own deep copy
+template <typename T>
+void NRMat<T>::copyonwrite()
+{
+#ifdef DEBUG
+        if (!count) laerror("Mat::copyonwrite of undefined matrix");
+#endif
+        if (*count > 1) {
+                (*count)--;
+                count = new int;
+                *count = 1;
+#ifdef MATPTR
+                T **newv = new T*[nn];
+                newv[0] = new T[mm*nn];
+                memcpy(newv[0], v[0], mm*nn*sizeof(T));
+                v = newv;
+                for (int i=1; i< nn; i++) v[i] = v[i-1] + mm;
+#else
+                T *newv = new T[mm*nn];
+                memcpy(newv, v, mm*nn*sizeof(T));
+                v = newv;
+#endif
+        }
+}
+
+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()");
+#endif
+        if (count)
+                if (*count > 1) {
+                        (*count)--;
+                        count = 0;
+                        v  = 0;
+                        nn = 0;
+                        mm = 0;
+                }
+        if (!count) {
+                count = new int;
+                *count = 1;
+                nn = n;
+                mm = m;
+#ifdef MATPTR
+                v = new T*[nn];
+                v[0] = new T[m*n];
+                for (int i=1; i< n; i++) v[i] = v[i-1] + m;
+#else
+                v = new T[m*n];
+#endif
+                return;
+        }
+        // At this point *count = 1, check if resize is necessary
+        if (n!=nn || m!=mm) {
+                nn = n;
+                mm = m;
+#ifdef MATPTR
+                delete[] (v[0]);
+#endif
+                delete[] v;
+#ifdef MATPTR
+                v = new T*[nn];
+                v[0] = new T[m*n];
+                for (int i=1; i< n; i++) v[i] = v[i-1] + m;
+#else
+                v = new T[m*n];
+#endif
+        }
+}
+
+
+
+
+
+
+
+
+
 // I/O
 template <typename T> extern ostream& operator<<(ostream &s, const NRMat<T> &x);
 template <typename T> extern istream& operator>>(istream  &s, NRMat<T> &x);

nonclass.h

@@ -3,7 +3,17 @@
 #include "mat.h"
 
 //MISC
-template <class T> extern const NRMat<T> diagonalmatrix(const NRVec<T> &x);
+export template <class T>
+const NRMat<T> diagonalmatrix(const NRVec<T> &x)
+{
+int n=x.size();
+NRMat<T> result((T)0,n,n);
+T *p = result[0];
+for(int j=0; j<n; j++) {*p  = x[j]; p+=(n+1);}
+return result;
+}
+
+//these just declared at the moment
 template <class T> extern const NRVec<T> lineof(const NRMat<T> &x, const int i); 
 template <class T> extern const NRVec<T> columnof(const NRMat<T> &x, const int i);
 template <class T> extern const NRVec<T> diagonalof(const NRMat<T> &x); 

smat.cc

@@ -33,59 +33,6 @@ NRSMat<T>::NRSMat(const NRMat<T> &rhs)
 }
 
 
-// dtor
-template <typename T>
-NRSMat<T>::~NRSMat()
-{
-	if (!count) return;
-	if (--(*count) <= 0) {
-		if (v) delete[] (v);
-		delete count;
-	}
-}
-
-
-// assignment with a physical copy
-template <typename T>
-NRSMat<T> & NRSMat<T>::operator|=(const NRSMat<T> &rhs)
-{
-	if (this != &rhs) {
-		if(!rhs.v) laerror("unallocated rhs in NRSMat operator |=");
-		if(count)
-			if(*count > 1) {	// detach from the other
-				--(*count);
-				nn = 0;
-				count = 0;
-				v = 0;
-			}
-		if (nn != rhs.nn) {
-			if(v) delete [] (v);
-			nn = rhs.nn;
-		}
-		if (!v) v = new T[NN2];
-		if (!count) count = new int;
-		*count = 1;
-		memcpy(v, rhs.v, NN2*sizeof(T));
-	}
-	return *this;
-}
-
-// assignment
-template <typename T>
-NRSMat<T> & NRSMat<T>::operator=(const NRSMat<T> & rhs)
-{
-	if (this == & rhs) return *this;
-	if (count)
-		if(--(*count) == 0) {
-			delete [] v;
-			delete count;
-		}
-	v = rhs.v;
-	nn = rhs.nn;
-	count = rhs.count;
-	if (count) (*count)++;
-	return *this;
-}
 
 // assing to diagonal
 template <typename T>
@@ -114,50 +61,7 @@ const T NRSMat<T>::trace() const
 	return tmp;
 }
 
-// make new instation of the Smat, deep copy
-template <typename T>
-void NRSMat<T>::copyonwrite()
-{
-#ifdef DEBUG
-	if (!count) laerror("probably an assignment to undefined Smat");
-#endif
-	if (*count > 1) {
-		(*count)--;
-		count = new int;
-		*count = 1;
-		T *newv = new T[NN2];
-		memcpy(newv, v, NN2*sizeof(T));
-		v = newv;
-	}
-}
 
-// resize Smat
-template <typename T>
-void NRSMat<T>::resize(const int n)
-{
-#ifdef DEBUG
-	if (n <= 0) laerror("illegal matrix dimension in resize of Smat");
-#endif
-	if (count)
-		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;
-		nn = n;
-		v = new T[NN2];
-		return;
-	}
-	if (n != nn) {
-		nn = n;
-		delete[] v;
-		v = new T[NN2];
-	}
-}
 
 // write matrix to the file with specific format
 template <typename T>

smat.h

@@ -315,6 +315,111 @@ inline NRSMat<T>:: operator const T*() const
 	return v;
 }
 
+//basic stuff to be available for any type ... must be in .h
+
+// dtor
+template <typename T>
+NRSMat<T>::~NRSMat()
+{
+        if (!count) return;
+        if (--(*count) <= 0) {
+                if (v) delete[] (v);
+                delete count;
+        }
+}
+
+// assignment with a physical copy
+template <typename T>
+NRSMat<T> & NRSMat<T>::operator|=(const NRSMat<T> &rhs)
+{
+        if (this != &rhs) {
+                if(!rhs.v) laerror("unallocated rhs in NRSMat operator |=");
+                if(count)
+                        if(*count > 1) {        // detach from the other
+                                --(*count);
+                                nn = 0;
+                                count = 0;
+                                v = 0;
+                        }
+                if (nn != rhs.nn) {
+                        if(v) delete [] (v);
+                        nn = rhs.nn;
+                }
+                if (!v) v = new T[NN2];
+                if (!count) count = new int;
+                *count = 1;
+                memcpy(v, rhs.v, NN2*sizeof(T));
+        }
+        return *this;
+}
+
+
+// assignment
+template <typename T>
+NRSMat<T> & NRSMat<T>::operator=(const NRSMat<T> & rhs)
+{
+        if (this == & rhs) return *this;
+        if (count)
+                if(--(*count) == 0) {
+                        delete [] v;
+                        delete count;
+                }
+        v = rhs.v;
+        nn = rhs.nn;
+        count = rhs.count;
+        if (count) (*count)++;
+        return *this;
+}
+
+
+// make new instation of the Smat, deep copy
+template <typename T>
+void NRSMat<T>::copyonwrite()
+{
+#ifdef DEBUG
+        if (!count) laerror("probably an assignment to undefined Smat");
+#endif
+        if (*count > 1) {
+                (*count)--;
+                count = new int;
+                *count = 1;
+                T *newv = new T[NN2];
+                memcpy(newv, v, NN2*sizeof(T));
+                v = newv;
+        }
+}
+
+
+// resize Smat
+template <typename T>
+void NRSMat<T>::resize(const int n)
+{
+#ifdef DEBUG
+        if (n <= 0) laerror("illegal matrix dimension in resize of Smat");
+#endif
+        if (count)
+                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;
+                nn = n;
+                v = new T[NN2];
+                return;
+        }
+        if (n != nn) {
+                nn = n;
+                delete[] v;
+                v = new T[NN2];
+        }
+}
+
+
+
 
 
 // I/O

sparsemat.h

@@ -197,11 +197,12 @@ r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)1); //saves a temporary and simplifi
 return r;
 }
 
+
 //add sparse to dense
 template<class T>
 NRMat<T> & NRMat<T>::operator+=(const SparseMat<T> &rhs)
 {
-if(nn!=rhs.nrows()||mm!=rhs.ncols()) laerror("incompatible matrices in +=");
+if((unsigned int)nn!=rhs.nrows()||(unsigned int)mm!=rhs.ncols()) laerror("incompatible matrices in +=");
 matel<T> *l=rhs.getlist();
 bool sym=rhs.issymmetric();
 while(l)
@@ -215,6 +216,7 @@ while(l)
 #endif
         l=l->next;
         }
+return *this;
 }
 
 

t.cc

@@ -8,7 +8,8 @@
 #include "fourindex.h"
 
 
-extern void test(const NRVec<double> &);
+extern void test(const NRVec<double> &x);
+
 
 double ad; 
 void f1(const double *c)
@@ -40,6 +41,8 @@ NRVec<double> x(1.,10);
 NRVec<double> y(2.,10);
 NRVec<double> z(-2.,10);
 
+if(0) test(x);
+
 y.axpy(3,x);
 
 y+=z;

t2.cc

@@ -0,0 +1,19 @@
+
+#include <time.h>
+#include "la.h"
+#include "traceback.h"
+#include "sparsemat.h"
+#include "matexp.h"
+#include "fourindex.h"
+
+
+void test(const NRVec<double> &x)
+{
+NRMat<double> aa(0.,2,2);
+NRMat<double> cc(aa);
+cc.copyonwrite(); cc[0][0]=1.;
+
+cout << aa << cc <<"\n";
+cout <<"test x= "<<x<<"\n";
+}
+

vec.cc

@@ -33,83 +33,8 @@ NRVec<T>::NRVec(const NRMat<T> &rhs)
 }
 #endif
 
-// dtor
-template <typename T>
-NRVec<T>::~NRVec()
-{
-	if(!count) return;
-	if(--(*count) <= 0) {
-		if(v) delete[] (v);
-		delete count;
-	}
-}
 
-// detach from a physical vector and make own copy
-template <typename T>
-void NRVec<T>::copyonwrite()
-{
-#ifdef DEBUG
-  if(!count) laerror("probably an assignment to undefined vector");
-#endif
-  if(*count > 1)
-  {
-    (*count)--;
-    count = new int;
-    *count = 1;
-    T *newv = new T[nn];
-    memcpy(newv, v, nn*sizeof(T));
-    v = newv;
-  }
-}
 
-// Asignment
-template <typename T>
-NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
-{
-  if (this != &rhs)
-  {
-    if(count)
-      if(--(*count) == 0)
-      {
-        delete[] v;
-        delete count;
-      }
-    v = rhs.v;
-    nn = rhs.nn;
-    count = rhs.count;
-    if(count) (*count)++;
-  }
-  return *this;
-}
-
-// Resize
-template <typename T>
-void NRVec<T>::resize(const int n)
-{
-#ifdef DEBUG
-  if(n<=0) laerror("illegal vector dimension");
-#endif
-  if(count)
-    if(*count > 1) {
-      (*count)--;
-      count = 0;
-      v = 0;
-      nn = 0;
-    }
-  if(!count) {
-    count = new int;
-    *count = 1;
-    nn = n;
-    v = new T[nn];
-    return;
-  }
-  // *count = 1 in this branch
-  if (n != nn) {
-    nn = n;
-    delete[] v;
-    v = new T[nn];
-  }
-}
 
 // ostream << NRVec
 template <typename T>
@@ -155,33 +80,6 @@ void NRVec<T>::fscanf(FILE *f, const char *format)
 			laerror("cannot read the vector eleemnt");
 }
 
-// assignmet with a physical copy
-template <typename T>
-NRVec<T> & NRVec<T>::operator|=(const NRVec<T> &rhs)
-{
-	if (this != &rhs) {
-#ifdef DEBUG
-		if (!rhs.v) laerror("unallocated rhs in NRVec operator |=");
-#endif
-		if (count)
-			if (*count > 1) {
-				--(*count);
-				nn = 0;
-				count = 0;
-				v = 0;
-			}
-		if (nn != rhs.nn) {
-			if (v) delete[] (v);
-			nn = rhs.nn;
-		}
-		if(!v) v = new T[nn];
-		if(!count) count = new int;
-		*count = 1;
-		memcpy(v, rhs.v, nn*sizeof(T));
-	}
-	return *this;
-}
-
 // unary minus
 template <typename T>
 const NRVec<T> NRVec<T>::operator-() const

vec.h

@@ -1,6 +1,8 @@
 #ifndef _LA_VEC_H_
 #define _LA_VEC_H_
 
+
+
 extern "C" {
 #include "cblas.h"
 }
@@ -44,6 +46,7 @@ template<class T> \
 #include "smat.h"
 #include "mat.h"
 
+
 // NRVec class
 template <typename T>
 class NRVec {
@@ -438,5 +441,113 @@ NRVECMAT_OPER2(Vec,-)
 
 // Few forward declarations
 
+//basic stuff which has to be in .h
+// dtor
+template <typename T>
+NRVec<T>::~NRVec()
+{
+        if(!count) return;
+        if(--(*count) <= 0) {
+                if(v) delete[] (v);
+                delete count;
+        }
+}
+
+// detach from a physical vector and make own copy
+template <typename T>
+void NRVec<T>::copyonwrite()
+{
+#ifdef DEBUG
+  if(!count) laerror("probably an assignment to undefined vector");
+#endif
+  if(*count > 1)
+  {
+    (*count)--;
+    count = new int;
+    *count = 1;
+    T *newv = new T[nn];
+    memcpy(newv, v, nn*sizeof(T));
+    v = newv;
+  }
+}
+
+// Asignment
+template <typename T>
+NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
+{
+  if (this != &rhs)
+  {
+    if(count)
+      if(--(*count) == 0)
+      {
+        delete[] v;
+        delete count;
+      }
+    v = rhs.v;
+    nn = rhs.nn;
+    count = rhs.count;
+    if(count) (*count)++;
+  }
+  return *this;
+}
+
+// Resize
+template <typename T>
+void NRVec<T>::resize(const int n)
+{
+#ifdef DEBUG
+  if(n<=0) laerror("illegal vector dimension");
+#endif
+  if(count)
+    if(*count > 1) {
+      (*count)--;
+      count = 0;
+      v = 0;
+      nn = 0;
+    }
+  if(!count) {
+    count = new int;
+    *count = 1;
+    nn = n;
+    v = new T[nn];
+    return;
+  }
+  // *count = 1 in this branch
+  if (n != nn) {
+    nn = n;
+    delete[] v;
+    v = new T[nn];
+  }
+}
+
+
+// assignmet with a physical copy
+template <typename T>
+NRVec<T> & NRVec<T>::operator|=(const NRVec<T> &rhs)
+{
+	if (this != &rhs) {
+#ifdef DEBUG
+		if (!rhs.v) laerror("unallocated rhs in NRVec operator |=");
+#endif
+		if (count)
+			if (*count > 1) {
+				--(*count);
+				nn = 0;
+				count = 0;
+				v = 0;
+			}
+		if (nn != rhs.nn) {
+			if (v) delete[] (v);
+			nn = rhs.nn;
+		}
+		if(!v) v = new T[nn];
+		if(!count) count = new int;
+		*count = 1;
+		memcpy(v, rhs.v, nn*sizeof(T));
+	}
+	return *this;
+}
+
+
 
 #endif /* _LA_VEC_H_ */