From e82adf91f898c664ec8b793dca712a01fb98591f Mon Sep 17 00:00:00 2001 From: jiri Date: Fri, 4 Sep 2009 08:09:32 +0000 Subject: [PATCH] *** empty log message *** --- bitvector.h | 2 +- mat.h | 15 +++++++++------ smat.h | 13 ++++++++----- sparsemat.h | 2 +- vec.h | 11 +++++++---- 5 files changed, 26 insertions(+), 17 deletions(-) diff --git a/bitvector.h b/bitvector.h index 4c4f368..d930d3b 100644 --- a/bitvector.h +++ b/bitvector.h @@ -53,7 +53,7 @@ public: void reset(const unsigned int i) {v[i/blockbits] &= ~(1<<(i%blockbits));}; const bool get(const unsigned int i) {return (v[i/blockbits] >>(i%blockbits))&1;}; const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;}; - void clear() {memset(v,0,nn*sizeof(bitvector_block));}; + void clear() {copyonwrite(); memset(v,0,nn*sizeof(bitvector_block));}; void fill() {memset(v,0xff,nn*sizeof(bitvector_block));}; bool operator!=(const bitvector &rhs) const; bool operator==(const bitvector &rhs) const {return !(*this != rhs);}; diff --git a/mat.h b/mat.h index 1fd2b11..a1ce209 100644 --- a/mat.h +++ b/mat.h @@ -42,6 +42,7 @@ public: inline NRMat(const T &a, const int n, const int m); NRMat(const T *a, const int n, const int m); inline NRMat(const NRMat &rhs); + NRMat(const typename LA_traits_complex::NRMat_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real explicit NRMat(const NRSMat &rhs); #ifdef MATPTR explicit NRMat(const NRVec &rhs, const int n, const int m, const int offset=0) :NRMat(&rhs[0][0] + offset ,n,m) {if (offset < 0 || n*m + offset > rhs.nn) laerror("matrix dimensions and offset incompatible with vector length");}; @@ -57,8 +58,7 @@ public: const bool operator==(const NRMat &rhs) const {return !(*this != rhs);}; inline int getcount() const {return count?*count:0;} NRMat & operator=(const NRMat &rhs); //assignment - void clear() {if(nn&&mm) LA_traits::clear((*this)[0],nn*mm);}; //zero out - void randomize(const T &x); //fill with random numbers + void randomize(const typename LA_traits::normtype &x); //fill with random numbers NRMat & operator=(const T &a); //assign a to diagonal NRMat & operator|=(const NRMat &rhs); //assignment to a new copy NRMat & operator+=(const T &a); //add diagonal @@ -88,6 +88,7 @@ public: const NRMat operator&(const NRMat &rhs) const; // direct sum const NRMat operator|(const NRMat &rhs) const; // direct product const NRVec operator*(const NRVec &rhs) const {NRVec result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec + const NRVec > operator*(const NRVec > &rhs) const {NRVec > result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec const NRVec rsum() const; //sum of rows const NRVec csum() const; //sum of columns const NRVec row(const int i, int l= -1) const; //row of, efficient @@ -95,6 +96,7 @@ public: const T* diagonalof(NRVec &, const bool divide=0, bool cache=false) const; //get diagonal void diagonalset(const NRVec &); //set diagonal elements void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const {r.gemv(beta,*this,trans,alpha,x);}; + void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const {r.gemv(beta,*this,trans,alpha,x);}; inline T* operator[](const int i); //subscripting: pointer to row i inline const T* operator[](const int i) const; inline T& operator()(const int i, const int j); // (i,j) subscripts @@ -105,6 +107,7 @@ public: void get(int fd, bool dimensions=1, bool transposed=false); void put(int fd, bool dimensions=1, bool transposed=false) const; void copyonwrite(); + void clear() {if(nn&&mm) {copyonwrite(); LA_traits::clear((*this)[0],nn*mm);}}; //zero out void resize(int n, int m); inline operator T*(); //get a pointer to the data inline operator const T*() const; @@ -118,7 +121,7 @@ public: const char transb, const T &alpha);//this = alpha*op( A )*op( B ) + beta*this void fprintf(FILE *f, const char *format, const int modulo) const; void fscanf(FILE *f, const char *format); - const double norm(const T scalar=(T)0) const; + const typename LA_traits::normtype norm(const T scalar=(T)0) const; void axpy(const T alpha, const NRMat &x); // this += a*x inline const T amax() const; const T trace() const; @@ -258,7 +261,7 @@ template inline T* NRMat::operator[](const int i) { #ifdef DEBUG - if (*count != 1) laerror("Mat lval use of [] with count > 1"); + if (_LA_count_check && *count != 1) laerror("Mat lval use of [] with count > 1"); if (i<0 || i>=nn) laerror("Mat [] out of range"); if (!v) laerror("[] for unallocated Mat"); #endif @@ -287,7 +290,7 @@ template inline T & NRMat::operator()(const int i, const int j) { #ifdef DEBUG - if (*count != 1) laerror("Mat lval use of (,) with count > 1"); + if (_LA_count_check && *count != 1) laerror("Mat lval use of (,) with count > 1"); if (i<0 || i>=nn &&nn>0 || j<0 || j>=mm && mm>0) laerror("Mat (,) out of range"); if (!v) laerror("(,) for unallocated Mat"); #endif @@ -626,7 +629,7 @@ public: inline T& operator() (const int i, const int j) { #ifdef DEBUG - if (*NRMat::count != 1) laerror("Mat lval use of (,) with count > 1"); + if (_LA_count_check && *NRMat::count != 1) laerror("Mat lval use of (,) with count > 1"); if (i<1 || i>NRMat::nn || j<1 || j>NRMat::mm) laerror("Mat (,) out of range"); if (!NRMat::v) laerror("(,) for unallocated Mat"); #endif diff --git a/smat.h b/smat.h index dfd7f91..4520303 100644 --- a/smat.h +++ b/smat.h @@ -38,12 +38,12 @@ public: inline NRSMat(const T &a, const int n); //Initialize to constant inline NRSMat(const T *a, const int n); // Initialize to array inline NRSMat(const NRSMat &rhs); // Copy constructor + NRSMat(const typename LA_traits_complex::NRSMat_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real explicit NRSMat(const NRMat &rhs); // symmetric part of general matrix explicit NRSMat(const NRVec &rhs, const int n); //construct matrix from vector NRSMat & operator|=(const NRSMat &rhs); //assignment to a new copy NRSMat & operator=(const NRSMat &rhs); //assignment - void clear() {LA_traits::clear(v,NN2);}; //zero out - void randomize(const T&x); + void randomize(const typename LA_traits::normtype &x); NRSMat & operator=(const T &a); //assign a to diagonal const bool operator!=(const NRSMat &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits::gencmp(v,rhs.v,NN2);} //memcmp for scalars else elementwise const bool operator==(const NRSMat &rhs) const {return !(*this != rhs);}; @@ -66,8 +66,10 @@ public: const T dot(const NRSMat &rhs) const; // Smat.Smat//@@@for complex do conjugate const T dot(const NRVec &rhs) const; //Smat(as vec).vec //@@@for complex do conjugate const NRVec operator*(const NRVec &rhs) const {NRVec result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec + const NRVec > operator*(const NRVec > &rhs) const {NRVec > result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec const T* diagonalof(NRVec &, const bool divide=0, bool cache=false) const; //get diagonal void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const {r.gemv(beta,*this,trans,alpha,x);}; + void gemv(const T beta, NRVec > &r, const char trans, const T alpha, const NRVec > &x) const {r.gemv(beta,*this,trans,alpha,x);}; inline const T& operator[](const int ij) const; inline T& operator[](const int ij); inline const T& operator()(const int i, const int j) const; @@ -76,13 +78,14 @@ public: inline int ncols() const; inline int size() const; inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant - const double norm(const T scalar=(T)0) const; + const typename LA_traits::normtype norm(const T scalar=(T)0) const; void axpy(const T alpha, const NRSMat &x); // this+= a*x inline const T amax() const; const T trace() const; void get(int fd, bool dimensions=1, bool transp=0); void put(int fd, bool dimensions=1, bool transp=0) const; void copyonwrite(); + void clear() {copyonwrite(); LA_traits::clear(v,NN2);}; //zero out void resize(const int n); inline operator T*(); //get a pointer to the data inline operator const T*() const; //get a pointer to the data @@ -283,7 +286,7 @@ template inline T & NRSMat::operator[](const int ij) { #ifdef DEBUG - if (*count != 1) laerror("lval [] with count > 1 in Smat"); + if (_LA_count_check && *count != 1) laerror("lval [] with count > 1 in Smat"); if (ij<0 || ij>=NN2) laerror("SMat [] out of range"); if (!v) laerror("[] for unallocated Smat"); #endif @@ -359,7 +362,7 @@ template inline T & NRSMat::operator()(const int i, const int j) { #ifdef DEBUG - if (*count != 1) laerror("lval (i,j) with count > 1 in Smat"); + if (_LA_count_check && *count != 1) laerror("lval (i,j) with count > 1 in Smat"); if (i<0 || i>=nn || j<0 || j>=nn) laerror("SMat (i,j) out of range"); if (!v) laerror("(i,j) for unallocated Smat"); #endif diff --git a/sparsemat.h b/sparsemat.h index be899f9..6633827 100644 --- a/sparsemat.h +++ b/sparsemat.h @@ -88,7 +88,6 @@ public: SparseMat & operator+=(const SparseMat &rhs); SparseMat & addtriangle(const SparseMat &rhs, const bool lower, const char sign); SparseMat & join(SparseMat &rhs); //more efficient +=, rhs will be emptied - void clear() {unsort();deletelist();} SparseMat & operator-=(const SparseMat &rhs); inline const SparseMat operator+(const T &rhs) const {return SparseMat(*this) += rhs;} inline const SparseMat operator-(const T &rhs) const {return SparseMat(*this) -= rhs;} @@ -133,6 +132,7 @@ public: inline bool issymmetric() const {return symmetric;} unsigned int length() const; void copyonwrite(); + void clear() {copyonwrite();unsort();deletelist();} void simplify(); const T trace() const; const T norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first diff --git a/vec.h b/vec.h index 4a61365..6ce250a 100644 --- a/vec.h +++ b/vec.h @@ -63,6 +63,7 @@ public: inline NRVec(const T *a, const int n); inline NRVec(T *a, const int n, bool skeleton); inline NRVec(const NRVec &rhs); + NRVec(const typename LA_traits_complex::NRVec_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real inline explicit NRVec(const NRSMat & S); #ifdef MATPTR explicit NRVec(const NRMat &rhs) : NRVec(&rhs[0][0],rhs.nrows()*rhs.ncols()) {}; @@ -71,8 +72,7 @@ public: #endif NRVec & operator=(const NRVec &rhs); NRVec & operator=(const T &a); //assign a to every element - void clear() {LA_traits::clear(v,nn);}; //zero out - void randomize(const T &x); + void randomize(const typename LA_traits::normtype &x); NRVec & operator|=(const NRVec &rhs); const bool operator!=(const NRVec &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits::gencmp(v,rhs.v,nn);} //memcmp for scalars else elementwise const bool operator==(const NRVec &rhs) const {return !(*this != rhs);}; @@ -99,6 +99,8 @@ public: void gemv(const T beta, const NRMat &a, const char trans, const T alpha, const NRVec &x); void gemv(const T beta, const NRSMat &a, const char trans /*just for compatibility*/, const T alpha, const NRVec &x); void gemv(const T beta, const SparseMat &a, const char trans, const T alpha, const NRVec &x,const bool treat_as_symmetric=false); + void gemv(const typename LA_traits_complex::Component_type beta, const typename LA_traits_complex::NRMat_Noncomplex_type &a, const char trans, const typename LA_traits_complex::Component_type alpha, const NRVec &x); + void gemv(const typename LA_traits_complex::Component_type beta, const typename LA_traits_complex::NRSMat_Noncomplex_type &a, const char trans, const typename LA_traits_complex::Component_type alpha, const NRVec &x); const NRVec operator*(const NRMat &mat) const {NRVec result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const NRSMat &mat) const {NRVec result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const SparseMat &mat) const {NRVec result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; @@ -115,11 +117,12 @@ public: void axpy(const T alpha, const NRVec &x); // this+= a*x void axpy(const T alpha, const T *x, const int stride=1); // this+= a*x void copyonwrite(); + void clear() {copyonwrite(); LA_traits::clear(v,nn);}; //zero out void resize(const int n); void get(int fd, bool dimensions=1, bool transp=0); void put(int fd, bool dimensions=1, bool transp=0) const; NRVec & normalize(); - inline const double norm() const; + inline const typename LA_traits::normtype norm() const; inline const T amax() const; inline const NRVec unitvector() const; void fprintf(FILE *f, const char *format, const int modulo) const; @@ -474,7 +477,7 @@ template inline T & NRVec::operator[](const int i) { #ifdef DEBUG - if(*count != 1) laerror("possible lval [] with count > 1"); + if(_LA_count_check && *count != 1) laerror("possible lval [] with count > 1"); if(i < 0 || i >= nn) laerror("NRVec out of range"); if(!v) laerror("[] on unallocated NRVec"); #endif