From badd89dcdb884fb4eb0ad38687e76c23747bf0bf Mon Sep 17 00:00:00 2001 From: jiri Date: Wed, 17 Mar 2004 16:39:07 +0000 Subject: [PATCH] *** empty log message *** --- mat.cc | 152 ----------------------------------------------- mat.h | 167 ++++++++++++++++++++++++++++++++++++++++++++++++++++ nonclass.h | 12 +++- smat.cc | 96 ------------------------------ smat.h | 105 +++++++++++++++++++++++++++++++++ sparsemat.h | 4 +- t.cc | 5 +- t2.cc | 19 ++++++ vec.cc | 102 -------------------------------- vec.h | 111 ++++++++++++++++++++++++++++++++++ 10 files changed, 420 insertions(+), 353 deletions(-) create mode 100644 t2.cc diff --git a/mat.cc b/mat.cc index b18c6d8..bafef93 100644 --- a/mat.cc +++ b/mat.cc @@ -14,43 +14,7 @@ template NRMat; * Templates first, specializations for BLAS next */ -// dtor -template -NRMat::~NRMat() -{ - if (!count) return; - if (--(*count) <= 0) { - if (v) { -#ifdef MATPTR - delete[] (v[0]); -#endif - delete[] v; - } - delete count; - } -} -// assign NRMat = NRMat -template -NRMat & NRMat::operator=(const NRMat &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 @@ -68,54 +32,7 @@ NRMat & NRMat::operator=(const T &a) return *this; } -// Explicit deep copy of NRmat -template -NRMat & NRMat::operator|=(const NRMat &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 @@ -215,76 +132,7 @@ const NRVec NRMat::rsum() const return result; } -// make detach Mat and make it's own deep copy -template -void NRMat::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 -void NRMat::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 diff --git a/mat.h b/mat.h index 2944f00..631e1aa 100644 --- a/mat.h +++ b/mat.h @@ -327,6 +327,173 @@ inline const complex NRMat< complex >::amax() const } +//basi stuff to be available for any type ... must be in .h +// dtor +template +NRMat::~NRMat() +{ + if (!count) return; + if (--(*count) <= 0) { + if (v) { +#ifdef MATPTR + delete[] (v[0]); +#endif + delete[] v; + } + delete count; + } +} + +// assign NRMat = NRMat +template +NRMat & NRMat::operator=(const NRMat &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 +NRMat & NRMat::operator|=(const NRMat &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 +void NRMat::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 +void NRMat::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 extern ostream& operator<<(ostream &s, const NRMat &x); template extern istream& operator>>(istream &s, NRMat &x); diff --git a/nonclass.h b/nonclass.h index de4b484..283c334 100644 --- a/nonclass.h +++ b/nonclass.h @@ -3,7 +3,17 @@ #include "mat.h" //MISC -template extern const NRMat diagonalmatrix(const NRVec &x); +export template +const NRMat diagonalmatrix(const NRVec &x) +{ +int n=x.size(); +NRMat result((T)0,n,n); +T *p = result[0]; +for(int j=0; j extern const NRVec lineof(const NRMat &x, const int i); template extern const NRVec columnof(const NRMat &x, const int i); template extern const NRVec diagonalof(const NRMat &x); diff --git a/smat.cc b/smat.cc index 9588ae8..b4ca579 100644 --- a/smat.cc +++ b/smat.cc @@ -33,59 +33,6 @@ NRSMat::NRSMat(const NRMat &rhs) } -// dtor -template -NRSMat::~NRSMat() -{ - if (!count) return; - if (--(*count) <= 0) { - if (v) delete[] (v); - delete count; - } -} - - -// assignment with a physical copy -template -NRSMat & NRSMat::operator|=(const NRSMat &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 -NRSMat & NRSMat::operator=(const NRSMat & 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 @@ -114,50 +61,7 @@ const T NRSMat::trace() const return tmp; } -// make new instation of the Smat, deep copy -template -void NRSMat::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 -void NRSMat::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 diff --git a/smat.h b/smat.h index c097128..0d860a0 100644 --- a/smat.h +++ b/smat.h @@ -315,6 +315,111 @@ inline NRSMat:: operator const T*() const return v; } +//basic stuff to be available for any type ... must be in .h + +// dtor +template +NRSMat::~NRSMat() +{ + if (!count) return; + if (--(*count) <= 0) { + if (v) delete[] (v); + delete count; + } +} + +// assignment with a physical copy +template +NRSMat & NRSMat::operator|=(const NRSMat &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 +NRSMat & NRSMat::operator=(const NRSMat & 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 +void NRSMat::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 +void NRSMat::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 diff --git a/sparsemat.h b/sparsemat.h index 96468ee..e05acb7 100644 --- a/sparsemat.h +++ b/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 NRMat & NRMat::operator+=(const SparseMat &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 *l=rhs.getlist(); bool sym=rhs.issymmetric(); while(l) @@ -215,6 +216,7 @@ while(l) #endif l=l->next; } +return *this; } diff --git a/t.cc b/t.cc index aca14ee..1e6d372 100644 --- a/t.cc +++ b/t.cc @@ -8,7 +8,8 @@ #include "fourindex.h" -extern void test(const NRVec &); +extern void test(const NRVec &x); + double ad; void f1(const double *c) @@ -40,6 +41,8 @@ NRVec x(1.,10); NRVec y(2.,10); NRVec z(-2.,10); +if(0) test(x); + y.axpy(3,x); y+=z; diff --git a/t2.cc b/t2.cc new file mode 100644 index 0000000..349f971 --- /dev/null +++ b/t2.cc @@ -0,0 +1,19 @@ + +#include +#include "la.h" +#include "traceback.h" +#include "sparsemat.h" +#include "matexp.h" +#include "fourindex.h" + + +void test(const NRVec &x) +{ +NRMat aa(0.,2,2); +NRMat cc(aa); +cc.copyonwrite(); cc[0][0]=1.; + +cout << aa << cc <<"\n"; +cout <<"test x= "<::NRVec(const NRMat &rhs) } #endif -// dtor -template -NRVec::~NRVec() -{ - if(!count) return; - if(--(*count) <= 0) { - if(v) delete[] (v); - delete count; - } -} -// detach from a physical vector and make own copy -template -void NRVec::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 -NRVec & NRVec::operator=(const NRVec &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 -void NRVec::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 @@ -155,33 +80,6 @@ void NRVec::fscanf(FILE *f, const char *format) laerror("cannot read the vector eleemnt"); } -// assignmet with a physical copy -template -NRVec & NRVec::operator|=(const NRVec &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 const NRVec NRVec::operator-() const diff --git a/vec.h b/vec.h index 33fb761..7b37d33 100644 --- a/vec.h +++ b/vec.h @@ -1,6 +1,8 @@ #ifndef _LA_VEC_H_ #define _LA_VEC_H_ + + extern "C" { #include "cblas.h" } @@ -44,6 +46,7 @@ template \ #include "smat.h" #include "mat.h" + // NRVec class template class NRVec { @@ -438,5 +441,113 @@ NRVECMAT_OPER2(Vec,-) // Few forward declarations +//basic stuff which has to be in .h +// dtor +template +NRVec::~NRVec() +{ + if(!count) return; + if(--(*count) <= 0) { + if(v) delete[] (v); + delete count; + } +} + +// detach from a physical vector and make own copy +template +void NRVec::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 +NRVec & NRVec::operator=(const NRVec &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 +void NRVec::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 +NRVec & NRVec::operator|=(const NRVec &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_ */